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I.  INTRODUCTION 

A.         BACKGROUND:  THE  MARITIME  SURVEILLANCE  ENVIRONMENT 

The  complexity  of  the  target  tracking  problem  varies  according  to  a  number  of 
factors.  Significant  research  has  been  devoted  to  the  study  of  trackers  in  varying  levels  of 
target  motion  model,  (plant)  noise  and  measurement  noise.  Performance  under  a  variety  of 
target  maneuvers  has  also  been  studied.  Clutter  environments  have  been  added  to  the 
model,  further  complicating  the  tracking  problem.  Tracking  problems  containing  multiple 
targets,  multiple  sensors,  and  target  maneuvers  are  most  complex,  and  extensive  research 
has  gone  into  developing  algorithms  for  these  situations.  A  maritime  surveillance 
application  for  such  a  tracker  is  the  focus  of  this  study.  The  characteristics  of  this 
environment  are  unique  in  many  respects,  presenting  an  interesting  and  challenging  tracking 
problem. 

One  of  the  most  problematic  aspects  of  maritime  surveillance  tracking  is  the  long 
time  interval  between  position  updates.  Typical  air  target  tracking  systems,  such  as  a  fire 
control  system,  will  receive  measurements  of  the  target  many  times  each  second  from 
tracking  radars  with  high  pulse  repetition  rates.  A  maritime  surveillance  tracking  system, 
however,  will  receive  measurements  of  contacts  at  intervals  measured  in  minutes  or  hours 
depending  on  such  factors  as  contact  geographic  location,  asset  availability,  and  source  level 
or  signal  strength.  Fortunately,  because  it  is  a  surveillance  system,  the  maritime 
surveillance  tracker  has  a  relaxed  accuracy  requirement  when  compared  to  an  air  target 


tracker:  accuracy  within  hundreds  or  a  few  thousand  yards  is  acceptable  in  this  application, 
but  clearly  unacceptable  in  air  target  tracking.  Also,  the  relatively  slow  speed  of  a  surface 
contact  mitigates  the  effect  of  longer  update  intervals.  Nonetheless,  the  undersampled 
environment  is  a  significant  challenge  that  a  maritime  surveillance  tracker  must  overcome. 

There  are  a  variety  of  participants  in  the  maritime  surveillance  effort:  active  radars 
from  ships  or  aircraft  at  sea,  nationally  tasked  sensors,  shipping  control  officials  using 
civilian  shipping  reports,  to  list  some.  While  the  sources  or  sensors  from  which  a  contact  s 
position  is  obtained  may  be  dissimilar  (acoustic  or  electronic  line  of  bearing,  active  radar,  or 
a  filed  plan  of  intended  movement),  the  data  is  integrated  into  the  maritime  surveillance 
picture  via  the  Sensor  Report.  From  this  report,  the  tracker  receives  either  a  line  of  bearing 
or  latitude  and  longitude  for  the  target  and  a  confidence  region.  This  confidence  region  may 
be  an  ellipse  (for  latitude  and  longitude  reports)  or  a  bearing  swath  (for  line  of  bearing 
reports)  and  is  sensor  dependent.  Other  information  that  may  be  contained  in  the  sensor 
report  includes  estimated  course  and  speed  and  feature  data  such  as  emitter  frequency  and 
pulse  repetition  frequency.  Sensor  reports  allow  dissimilar  sensors  to  participate  in  the 
tracking  effort. 

Unlike  many  tracking  problems  where  position  information  is  received  at  regular 
intervals,  there  is  no  synchronization  of  sensor  reports  in  maritime  surveillance.  The 
various  sensors  participating  each  operate  according  to  their  own  reporting  abilities  and 
frequently  have  irregular  reporting  periods.  Consequently,  batch  processing  of 
measurements  is  impractical. 


Finally,  the  number  of  tracks  that  a  maritime  surveillance  tracker  must  handle  to  be 
of  value  further  complicates  the  design.  With  the  goal  of  providing  situational  awareness  of 
the  surface  traffic  in  a  geographic  area  to  naval  forces,  a  maritime  surveillance  tracker 
would  have  to  support  thousands  of  tracks.  While  advances  in  computing  and  memory 
storage  technologies  will  enable  the  design  of  a  high  volume  tracker,  the  number  of  contacts 
to  be  tracked  in  a  geographic  area  will  vary  with  time  such  that  the  tracking  system  is  faced 
with  an  unknown  number  of  targets.  This  restriction  eliminates  many  tracking  algorithms 
from  consideration. 


B.         THE  TRACK  SPLIT  APPROACH 

The  track  split  approach  is  a  multi-target  tracking  algorithm  that  lends  itself  to 
several  problems,  including  single  target  in  clutter,  multiple  targets  in  clutter,  and  multiple 
targets  in  clutter  with  the  number  of  targets  unknown.  Unlike  some  simple  trackers  which 
assume  the  measurement  closest  to  the  estimate  is  the  correct  one,  track  splitting  trackers 
will  use  all  measurements  that  fall  within  a  validation  region,  or  gate,  to  update  the  target 
state.  A  likelihood  function  is  then  used  to  eliminate  unlikely  sequences  that  presumably 
represent  false  tracks.  In  this  way,  the  tracker  postpones  important  decisions  about  which 
measurements  originated  from  a  target  until  additional  measurement  information  is 
received. 

The  track  split  approach  shares  the  Kalman  Filter  based  prediction  and  update 
sequence  used  by  other  trackers.  This  widely  used  technique  obtains  the  optimal  estimate  of 
the  target  state,  for  a  given  dynamic  model,  by  recursively  minimizing  the  mean  square 


error.  The  filter  is  tunable  to  a  variety  of  applications  and  environments  by  properly 
selecting  the  dynamic  and  noise  model.  Asynchronous  or  missed  measurements  are  readily 
handled  by  the  filter  and  a  measure  of  tracker  accuracy  (relative  to  the  specified  model)  is 
always  available  in  the  form  of  the  target's  covariance  matrix.  Because  of  its  versatility,  a 
Kalman  Filter  based  Track  Splitting  tracker  can  be  developed  for  a  maritime  surveillance 
role. 

The  track  split  approach  was  selected  for  this  study  because  it  is  the  simplest 
algorithm  capable  of  meeting  the  requirements  imposed  by  the  maritime  surveillance 
environment.  Track  splitting  is  a  technique  for  allowing  the  Kalman  Filter  to  handle  false 
measurements  and  multiple  contacts.  Furthermore,  this  approach  does  not  require  the 
number  of  targets  to  be  known  and  will  automatically  initialize  new  tracks  when  more  than 
one  measurement  falls  within  the  gate  of  a  currently  held  track.  Although  other  algorithms 
are  also  capable  of  performing  under  these  conditions  (most  notably,  the  Multiple 
Hypothesis  Tracker  (MHT)),  the  track  split  approach  is  a  logical  starting  point  in  the  design 
and  evaluation  of  a  maritime  surveillance  tracker. 


C.         THESIS  OUTLINE 

The  remainder  of  this  thesis  is  organized  as  follows:  Chapter  II  will  discuss  the 
contact  and  measurement  model  by  presenting  the  state  space  dynamic  representation  and 
noise  model.  Chapter  IH  will  detail  the  tracker  algorithm,  including  a  review  of  standard 
and  extended  Kalman  Filter  equations,  track  split  logic,  gating,  and  estimate  error  ellipses. 
Chapter  IV  will  present  and  discuss  the  simulation  parameters  and  results.   Chapter  V  will 
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draw  conclusions  about  this  study  and  outline  some  areas  for  future  research.  The  tracker 
simulation  code,  written  for  MATLAB®,  is  contained  in  the  Appendix. 


II.    CONTACT  AND  MEASUREMENT  MODEL 


A. 


CONTACT  MOTION  MODEL 


The  tracker  assumes  contacts  have  a  two-dimensional  (x-y),  second  order  (position 
and  velocity)  dynamic  model.  The  model  is  implemented  in  discrete  form  according  to: 


x  t  +  1    =    F*   •  x  k    +    v  k 


2.1 


where  x  is  the  state  space  vector  composed  of  x  position,  x  velocity,  y  position,  and  y 
velocity.  F  is  the  linear  matrix: 


F*  = 


1 
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0 
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Here,  AT  is  the  time  interval  between  samples.  Plant  noise  is  included  in  the 
dynamic  model  to  account  for  contact  deviations  from  straight-line  motion.  Uk  is  a  zero 
mean  Gaussian  random  variable  with  known  covariance,  Q  ,  given  by: 
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2.3 


The  parameter  q"  is  used  to  adjust  the  magnitude  of  plant  noise.  Properly  selecting 
this  parameter  allows  the  designer  to  optimize  tracker  performance  for  the  expected 
behavior  of  the  contacts.  Contacts  that  maintain  relatively  straight  line  motion,  in  open 
ocean  transit  for  example,  will  not  deviate  greatly  from  the  dynamic  model.  A  low  value 
for  q2  will  give  optimum  tracker  performance.  Maneuvering  contacts,  however,  conform  to 
the  model  poorly.  A  large  value  for  q2  allows  the  tracker  to  maintain  track. 


B.         MEASUREMENT  MODEL 

Sensors  feeding  into  the  maritime  surveillance  tracking  system  obtain  contact 
positions  as  either  single  dimensional  lines  of  bearing  or  two  dimensional  range  and  bearing 
measurements.  Passive  systems,  such  as  acoustic  hydrophones  or  electronic  direction 
finding  receivers,  will  produce  lines  of  bearing.  Active  radars  and  electronic  geolocation 
systems  employing  techniques  such  as  Time  Difference  of  Arrival  measurements  produce 
two-dimensional  position  reports.  All  sensors  have  an  associated  measurement  error  which 
is  modeled  as  a  zero  mean  Gaussian  random  variable. 

Position  reports  contain  the  sensor's  measurement,  along  with  a  metric  of 
measurement  uncertainty.  In  the  case  of  Line  of  Bearing  reports,  a  95%  confidence  swath 
represents  about  2.5  standard  deviations  of  the  measurement  noise.  Thus,  Gq  can  be  readily 
found  from  the  sensor  report.  Latitude  and  Longitude  reports  will  report  measurement  error 
statistics  in  the  form  of  an  error  ellipse.  This  ellipse  is  related  to  the  eigenvalues  and 
eigenvectors  of  the  sensor's  covariance  matrix,  R,  (Kirk,  1975): 
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where 

X\  is  the  eigenvalue  associated  with  the  X'  eigenvector  of  R 
Xi  is  the  eigenvalue  associated  with  the  Y'  eigenvector  of  R 


c2  is  a  constant  that  determines  the  size  of  the  ellipse 


The  relationship  between  c  and  the  probability  of  a  measurement  lying  within  the 
ellipse  is  tabulated  for  the  two  dimensional  case  (Kirk,  1975).  In  this  study,  c  is  chosen  to 
be  6,  representing  approximately  95%  probability  that  the  measurement  falls  within  the 
error  ellipse. 

The  tracker  must  invert  the  ellipse  generating  process  to  obtain  the  measurement 

covariance  matrix.  The  sensor  report  contains  the  semi-major  axis  length  (a),  semi-minor 

axis  length  (b),  and  orientation  of  the  ellipse  ($).  The  eigenvalues  can  be  found  by: 


K 

= 

a 
2c 

K 

— 

b 

2.5 

2c 

The  orientation  of  the  semi-major  and  semi-minor  axes  may  be  expressed  as  unit 
vectors,  giving  the  eigenvectors  of  the  matrix  that  generated  the  ellipse  directly.    A  similar 


matrix  to  R  (one  with  the  same  characteristic  equation)  can  be  obtained  as  follows  (Kaplan, 
1981): 


R'  =  C  *  B  *  C  A 

where  C  is  composed  of  the  eigenvectors  obtained  from  the  ellipse  orientation: 

C  =  [X'    Y'] 


2.6 


2.7 


and  B  is  composed  of  the  eigenvectors  obtained  from  the  semi-major  and  semi-minor  axes: 


B  = 


A,      0 

0      I, 


2.8 


While  the  matrix  R'  is  not  necessarily  identical  to  the  sensor  s  covariance  matrix,  it  has  the 
same  characteristic  equation  and  is  therefore  acceptable  for  use  in  the  tracker. 

A  coordinate  transformation  from  polar  (range  and  bearing)  to  cartesian  (latitude 
and  longitude)  occurs  at  the  sensor.  For  purposes  of  this  simulation,  the  cartesian  position 
in  distance  units  is  transformed  into  latitude  and  longitude  coordinates  using  a  flat  earth 
model,  with  one  degree  equal  to  60  nautical  miles  (nm).  A  nonlinear  rotation  is  used  to 
obtain  the  measurement  in  cartesian: 


X 
Y 


r„  sin 
r_  cos 


2.9 


where  rm  and  6m  are  the  measured  range  and  bearing  to  the  contact.  Qm  is  the  angle  measured 
from  a  vertical  axes,  or  "True  North"  reference.  The  measurement  covariance  matrix  must 
also  be  transformed  from  polar  to  cartesian  coordinates.  The  resulting  elements  of  the 
cartesian  covariance  matrix  (Bar  Shalom  and  Li,  1995): 
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where  2 10 

*n  =  r>.2sin  26m  +  <r;cos  2em 

Rn  =  r^l  cos  2  #„  +  *  ?  sin  2  #m 
^12  =  f^r  "  C^^sin  #m  cos  #m 


where  oT  and  <7<9  are  the  range  and  bearing  standard  deviations  of  the  sensor.  For  the 
simulation,  the  measurement  and  covariance  are  transformed  from  polar  to  cartesian 
coordinates  by  the  m-file  pole2cart.  The  listing  is  contained  in  the  Appendix. 

Equipped  with  these  models,  a  tracker  can  be  implemented  that  will  produce 
estimates  from  line  of  bearing  or  cartesian  position  reports.  The  measurement  and  sensor 
covariance  can  be  determined  either  directly,  as  in  the  case  of  the  line  of  bearing  report,  or 
by  analyzing  the  error  ellipse  associated  with  a  cartesian  report.  Since  measurements  are 
generated  for  the  simulation  in  this  thesis,  a  technique  for  linearizing  range  and  bearing 
measurements  to  obtain  the  position  in  cartesian  coordinates  is  also  required.  The 
development  of  a  track  splitting  tracker  in  the  upcoming  chapter  will  draw  on  the  models 
discussed  here. 
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III.       THEORETICAL  BASIS  FOR  THE  TRACK  SPLITTING  TRACKER 

A.         THE  KALMAN  FILTER 

The  Kalman  Filter  is  a  recursive  algorithm  for  estimating  the  contact's  current  state 
vector  from  its  past  estimate  and  current  noisy  measurement.  With  accurately  modeled 
plant  dynamics  and  sensor  noise,  the  algorithm  minimizes  the  mean  square  error  in  the 
resulting  estimates.  Several  excellent  sources  exist  for  detailed  derivations  of  the  Kalman 
Filter  (Bar  Shalom  and  Li,  1998).  Only  the  resulting  equations  are  presented  here.  The 
following  definitions  and  notations  will  be  used: 

xfc+1,j,  =  the  estimate  of  state  vector  (x)  at  time  k  + 1  given  the  measurement 

at  time  k  (prediction) 
xklk    =  the  estimate  of  state  vector  (x)  at  time  k  given  the  measurement 

at  time  k  (corrected  estimate) 
P*+n*  - tne  covariance  of  state  vector  (x)  at  time  k  + 1  given  the  measurement  at  time  k 
Fklk    =  the  covariance  of  state  vector  (x)  at  time  k  given  the  measurement  at  time  k 

The  Kalman  Filter  recursion  contains  a  prediction  and  a  correction  step.  The 
mechanics  of  these  steps  are  now  described. 

Prediction:  The  weighted  least  squares  estimate  and  associated  covariance  matrix 
for  the  next  sample  time  are  predicted  based  on  plant  and  noise  models  according  to  these 
equations: 

X   k  +  1  \k      =      *  k  X   Jt|Jt 

/  D.  1 

P  =   F   P      F      +  O 

*   k  +  ilk  r  kr  k\k   r  k     T    v   k 
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3.2 


Correction:  The  estimate  and  covariance  matrix  at  time  k+1  are  corrected  with  the 
new  (time  k+1)  measurement. 

P*+il»+i  =  (I  ~  KA+,H)Pt+1|Jfc  (I  -  Ki+1H)'  +  K^R^K^ 

where 

Z*+l  —  ZA+1  _  "*k+\\k 

K*+i  =  P*+n*HTHPA+1|tH'+RJfc+ir1 

The  standard  Kalman  Filter  presented  above  can  be  used  when  the  plant  and 
measurement  models  are  linear,  as  in  the  case  of  cartesian  position  reports.  The 
measurement  matrix,  H,  for  the  standard  filter  is: 


H  = 


10    0    0" 
0    0     10 


3.3 


The  line  of  bearing  report,  however,  requires  a  nonlinear  measurement  matrix  since 


9(x)     =  arc  tan 


'O 


KX*J 


3.4 


The  Extended  Kalman  Filter  provides  a  means  to  handle  nonlinear  state  or  measurement 
equations  by  linearizing  the  nonlinear  equations  about  the  estimate.  For  the  line  of  bearing 
problem,  the  linearized  measurement  matrix  obtained  from  a  first  order  Taylor  Series 


expansion  is: 

H,  = 

xi 

0 

-    *3 

2     ,          2 

2               2 

*-**+i|* 


3.5 
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The  measurement  matrix  is  linearized  about  the  predicted  state  and  must  be  calculated  with 
each  new  measurement.  Once  the  linearized  measurement  matrix  is  obtained,  it  may  be 
used  in  the  correction  step  in  the  same  manner  that  the  static  measurement  matrix  (Eqn.  3.3) 
is  used  with  the  standard  Kalman  Filter  to  obtain  the  corrected  estimate. 

Since  the  Kalman  Filter  is  recursive,  an  initial  state  estimate  and  state  covariance  are 
required  to  begin  the  prediction  and  correction  process.  A  first  order  polynomial  curve  fit 
method  for  obtaining  these  from  two  (cartesian)  measurements  may  be  used  (Bar  Shalom 
and  Li,  1998).  The  latest  measurement  is  used  for  the  position  estimate,  while  the  velocity 
estimate  is  obtained  from  the  difference  in  position.  The  initial  covariance  is  calculated: 

l 


Pn  = 


1 

0 

0 

0 

1 

- 1 

0 

0 

AT 

AT 

0 

1 

0 

0 

1 

-  1 

0 

0 

AT 


AT 


R 


cart 

0 


0 


R 


cart 


1 

AT 

0  0 

-  1 


0 


AT 
0  0 


0 

0 

1 

1 

AT 

0 

0 

0 

-  1 

AT  . 
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Rcart  is  as  derived  in  the  previous  chapter  (Eqn.  2.10),  with  the  subscript  denoting  the 
measurement  with  which  the  covariance  matrix  is  associated.  The  m-file  kaI2init  generates 
the  filter  initializations  from  the  plant  model  and  simulation  parameters. 

The  Kalman  Filter  recursions  form  the  backbone  of  the  tracker  algorithm,  predicting 
and  correcting  the  contact's  state  estimate  based  on  the  most  recent  state  estimate,  the  plant 
model,  and  the  most  recent  measurement  of  the  contact.  In  the  multitarget  tracking 
problem,  some  data  association  decisions  must  be  made  to  pair  measurements  with  existing 
tracks  before  this  process  can  begin.  A  technique  for  comparing  the  features  of  the 
measurement  with  the  known  features  of  the  target  can  be  employed  to  quickly  reduce  the 
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number  of  candidates  for  update.  Also,  a  gating  procedure  will  eliminate  the  most  unlikely 
candidates  by  considering  only  the  measurements  within  a  region  around  the  prediction. 


B.         MEASUREMENT  REJECTION  USING  FEATURES 

One  technique  employed  in  this  study  relies  on  the  assumption  that  the  features  for 
both  a  contact  being  tracked  and  a  false  measurement  are  normally  distributed  with  different 
statistics.  The  vector  of  features  used  here  is: 


y  = 


emitter  frequency  (GHz) 
emitter  pulse  repetition  frequency  (pps) 


3.7 


Clearly,  not  all  sensors  would  include  these  features  in  the  position  report,  so  clutter 
rejection  using  features  can  not  always  be  accomplished.  The  statistics  of  the  feature 
random  variables  are  denoted: 


y  = 


E[frequency 
E[prf] 


3.8 


frequency 


Distances  are  then  formed  to  compare  each  measurement  with  statistics  of  the  contact  and 
clutter  models.    Here,  the  subscript  "7"  denotes  target,  and  "C"  denotes  clutter,  or  false 


target: 


d  T  =  (y  -  y  r  /£  rVy  -  y  T ) 
dc  =  (y  -  Tc  )'i~cl(J  -  Tc  ) 


3.9 
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Small  distances  indicate  a  measurement  that  closely  matches  the  feature  model,  based  on 
the  feature  statistics.  Forming  a  ratio  and  establishing  a  threshold  provides  a  criterion  for 
rejecting  a  measurement  as  clutter: 

■^-Z—  >    Threshold  3.10 

d  c 

Alternatively,  when  the  clutter  feature  means  are  not  known,  a  gating  technique  can 
be  used  with  the  target  feature  distance  alone.  In  this  case,  measurements  are  rejected  when 
dr  exceeds  a  Chi-Squared  distribution  threshold.  Using  a  threshold  of  6  will  retain  the 
correct  measurement  with  95%  confidence  while  discarding  measurements  with  features 
that  conform  poorly  with  the  expected  values.  With  both  feature  rejection  techniques, 
measurements  that  meet  the  feature  criteria  remain  valid  candidates  for  update.  Further 
processing  will  reduce  the  number  of  candidates  based  on  the  state  estimate  and  position 
measurement. 


C.        ELLIPSOIDAL  GATING 

Just  as  a  confidence  ellipse  exists  around  a  noisy  measurement,  a  region  can  be 
defined  around  an  estimate  that  will  establish  the  area  within  which  a  candidate  observation 
must  lie  with  some  probability.  The  smallest  region  for  a  given  probability  is  an  ellipse 
defined  by  the  innovation  and  innovation  covariance  resulting  from  a  measurement,  (Bar 
Shalom  and  Li,  1995): 
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t-^nJ^'fc-^nJ^y 


where  3.11 

L  =  HP;|i  ,H'+  R  =  the  innovation    covariance 

Y  =  chi  -  squared    threshold 

z  -  ik+lVc  =  the  measuremen    t's  innovation 

The  above  quadratic  (norm  of  the  innovation  squared)  is  a  Chi-Squared  random  variable 
with  the  number  of  degrees  of  freedom  equal  to  the  dimension  of  the  measurement  (two  for 
cartesian,  one  for  line  of  bearing).  Thus,  for  each  measurement  received  and  each  track 
held,  the  measurement  innovation  and  innovation  covariance  must  be  calculated.  With 
these,  the  norm  of  the  innovation  squared  may  be  calculated  and  compared  to  the  threshold. 
Norms  exceeding  the  threshold  indicate  a  pairing  that  is  highly  unlikely  and  need  not  be 
considered. 

The  application  of  this  gating  quadratic  is  straightforward  for  the  cartesian  position 
report.  In  these  simulations,  a  threshold  value  of  6  has  been  used  for  cartesian  reports, 
representing  a  95%  confidence  gate.  For  line  of  bearing  reports,  the  linearized  measurement 
matrix  must  be  used  in  determining  the  innovation  covariance.  The  95%  threshold  for  this 
single  degree  of  freedom  distribution  has  been  approximated  as  4  in  this  simulation,  but  in 
reality  is  slightly  less. 


D.         TRACK  SPLITTING  AND  THE  LIKELIHOOD  FUNCTION 

If  only  one  measurement  remains  a  candidate  for  track  updating  after  processing 
with  feature  measurements  and  gates,  no  further  tracker  logic  is  required.    The  candidate 
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measurement  is  passed  off  to  the  Kalman  Filter  (standard  or  extended)  and  the  corrected 
estimate  can  be  obtained.  If  no  measurements  remain  a  candidate  after  processing,  options 
can  be  included  in  the  tracker  to  increase  thresholds  in  the  feature  or  ellipsoidal  gate 
algorithms  in  an  effort  to  obtain  a  candidate.  (No  such  logic  has  been  included  in  this 
simulation,  however).  When  more  than  one  measurement  remains  a  candidate,  the  tracker 
is  forced  to  make  a  decision  as  to  which  one  most  likely  originated  from  the  contact  being 
tracked.  With  the  track  split  approach,  this  decision  is  deferred  until  additional  information 
(measurements)  can  be  obtained.  The  predicted  state  is  updated  with  each  measurement  (by 
the  Kalman  Filter  correction  step)  so  that  a  single  track  will  become  n  tracks  after 
processing,  where  n  is  the  number  of  candidate  measurements.  For  each  of  these  tracks,  a 
Kalman  Filter  prediction  will  be  obtained  so  that  a  gate  may  be  established  for  subsequent 
measurements.  Feature  processing,  gating,  splitting,  and  updating  will  be  repeated  for  each 
track  as  new  measurements  are  received. 

With  this  logic,  the  number  of  tracks  held  in  the  system  can  grow  exponentially.  To 
aid  in  identifying  the  correct  track,  (that  is,  the  one  most  closely  matching  the  plant  model)  a 
likelihood  function  based  on  the  Gaussian  assumptions  of  the  plant  and  measurement  model 
is  computed  recursively  each  time  a  measurement  updates  a  track  (Bar  Shalom  and  Li, 
1998): 

^=VI  +  [z-zi+j2-1[z-zi+1J  3.12 

That  is,  the  likelihood  that  a  sequence  of  measurements  originated  from  a  track  conforming 
to  the  model  is  the  sum  of  the  previous  score  and  the  norm  of  the  innovation  squared.   A 
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low  score  indicates  that  the  sequence  conforms  to  the  model  well  and  is  likely  a  valid  track. 
Thus,  a  metric  exists  to  aid  in  managing  the  number  of  tracks  in  the  system 


E.         TRACK  DROPPING  CRITERIA 

A  Track  Dropping  algorithm  is  essential  for  a  track  splitting  tracker  to  prevent 
overloading  of  the  computational  and  data  storage  capacities  of  the  system.  Track  scores, 
computed  by  the  likelihood  function  (Eqn.  3.12),  are  one  available  tool  for  selecting  the 
state  estimates  that  have  arisen  out  of  the  most  unlikely  measurement  sequences.  By 
choosing  an  appropriate  threshold,  the  number  of  false  tracks  held  in  the  system  can  be  kept 
to  a  minimum  without  inadvertantly  dropping  valid  tracks.  Another  parameter  that  aids  in 
track  management  is  the  most  recent  time  at  which  an  estimate  was  updated  with  a 
measurement,  or  update  time.  A  threshold  can  be  set  for  this  parameter  as  well,  triggering 
the  tracker  to  automatically  eliminate  a  held  track  once  it  is  exceeded.  A  combination  of 
these  techniques  is  included  in  this  simulation. 
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IV.   COMPUTER  SIMULATION 

A.    OVERVIEW 

The  listing  for  a  simulation  of  a  track  splitting  tracker  written  for  MATLAB®  is 
contained  in  the  Appendix.  This  tracker  incorporates  the  theoretical  models  previously 
discussed.  Figure  1  shows  the  functional  flow  of  the  tracker  simulation.  Two  dimensional 
constant  velocity  motion  is  simulated  for  two  targets  in  a  flat  earth  coordinate  system.  An 
instantaneous  turn  may  be  specified  by  selecting  a  time  parameter  for  that  occurrence.  Two 
attributes,  emitter  frequency  and  pulse  repetition  frequency  (PRF),  may  also  be  associated 
with  the  contacts.  A  range  and  bearing  measurement  of  contact  position  from  a  fixed  sensor 
position  is  made  nearly  every  90  minutes.  The  contact's  features  are  also  obtained  at  each 
measurement.  Bearing-only  measurement  times  may  also  be  specified.  Along  with  the 
noisy  measurements  of  contact  position,  false  measurements  that  did  not  originate  from  the 
contact  are  generated  at  sample  times  in  an  area  around  the  contact's  position.  These  false 
measurements  have  their  own  associated  feature  measurements.  Thus,  a  position 
measurement  and  the  associated  feature  measurements  at  a  specific  measurement  time  can 
be  considered  a  single  sensor  report. 

After  generating  all  measurements  at  a  sample  time,  track  processing  begins.  The 
feature  rejection  loop  compares  the  expected  feature  values  for  a  given  track  with  the 
measured  values  for  each  report.  Reports  that  do  not  conform  well  to  the  track's  expected 
features  are  removed  from  the  set  of  candidates  for  the  track  being  processed,  but  will  be 
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Figure  1.  Functional  Flow  of  Track-Splitting  Simulation 
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considered  for  subsequent  tracks.  The  remaining  candidate  reports  are  filtered  through  a 
geometric  gate  around  the  track's  expected  position.  For  all  measurements  that  fall  within 
the  gate,  a  corrected  state  estimate  and  likelihood  score  are  computed.  If  no  measurement 
falls  within  the  gate,  the  track  is  not  updated  at  that  sample  time.  Once  all  tracks  have  been 
processed  with  the  measurements  for  a  given  sample  time,  a  track  management  routine 
begins.  This  routine  performs  two  functions:  track  dropping  and  track  association.  The 
track  dropping  function  checks  likelihood  scores  and  estimate  update  times.  If  these  exceed 
thresholds,  the  tracks  are  dropped  and  will  not  be  considered  in  future  processing.  The  track 
association  function  compares  the  scores  of  all  estimates  that  are  candidates  for  the  original 
two  targets,  and  selects  the  estimates  that  correspond  with  the  lowest  scores. 

After  looping  through  all  measurement  times,  the  simulation  generates  several 
output  plots.  True  contact  motion,  noisy  measurements,  false  measurements,  predicted 
position  estimates,  and  corrected  position  estimates  are  displayed  in  several  figures.  Also,  a 
statistics  matrix  that  is  useful  in  tracker  performance  evaluation  is  displayed  in  the 
MATLAB  command  window.  Items  such  as  the  number  of  tracks  held,  the  number  of 
clutter  points  generated,  and  the  measurement  and  estimate  error  at  each  sample  time  can  be 
reviewed. 


B.         SIMULATION  1:  FEATURE  REJECTION  TEST 

The  first  simulation  demonstrates  the  ability  of  the  feature  rejection  algorithm 
discussed  in  Section  B  of  Chapter  III  to  discriminate  between  false  measurements  and  target 
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measurements.  Here,  clutter  frequency  and  PRF  means  are  well  separated  from  the  two 
targets  relative  to  the  standard  deviation  of  the  sensor's  measurements.  The  parameters  for 
the  simulation  are  as  follows: 

Simulation: 


number  of  measurements :  1 6 
clutter  density:  5  false  targets  per  60nm  X  60nm  area 
Sensor: 

Position:  18.00N         113.00W 

ur:  5  nm 

<Je-  0.5  degrees 

Gfreq.    0.001  GHZ 
Gprf.    0.1  ppS 

Mean  Clutter  Features: 

f clutter-     1.4GHZ 
Prfdutter-    3.8  ppS 

Target  1: 

Initial  Position:  22.00N  110.00W 

Heading:  145T 
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Speed:  13  Kts 

Turns:  To  180T  at  590  minutes  in  simulation 

fmean'.    1.406  GHz 

Prfmean-    4.4  pps 

q2:  lxlO"8 
Target  2: 

Initial  Position:  22.00N  109.00W 

Heading:  120T 
Speed:  15  Kts 
Turns:  none 

fmean-.    1.394  GHz 

prfmean-    3.2  pps 

q2:   lxlO"8 

Measurements:  14  cartesian  measurements  at  about  90  minute  intervals.  Additionally,  one 
Line  of  Bearing  measurement  will  be  obtained  30  minutes  before  the  5th  cartesian 
measurement.  A  second  line  of  bearing  measurement  will  be  obtained  45  minutes 
before  the  10th  cartesian  measurement. 
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Filter  Settings  and  Thresholds: 

Geographic  gating:  y  =  6,  corresponding  to  95%  probability  for  a  chi-squared 
distributed  random  variable. 

Delete  track  score  threshold:  48 
Delete  track  time-late  threshold:  200  minutes 

Figure  2  shows  the  true  track  for  targets  one  and  two.  Target  one  will  follow  the 
course  change  model,  and  target  two  will  maintain  a  straight  line  track. 
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Figure  2.  True  Target  Motion  for  Targets  1  and  2. 
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Figure  3  shows  the  true  target  position  at  sampling  times,  target  and  false 
measurements,  and  tracker  position  predictions  (z^.;).  The  ellipses  represent  the  95% 
confidence  regions  around  the  predictions  and  are  obtained  from  J*k\k-i  ■  These  ellipses  are 
related  to  the  geographic  gate  used  by  the  tracker.  It  is  noted  that  many  measurements  are 
geographically  feasible  candidates,  but  it  will  be  shown  that,  for  this  simulation,  the  feature 
rejection  loop  effectively  removes  most  from  consideration. 

Target  Motion  and  Measurements  at  Sample  Times 


-113     -112      -111      -110      -109      -108      -107      -106      -105      -104 

Deg  Lon 

Figure  3.  True  Target  Position  and  Measurements  At  Sample  Times. 
The  ellipses  are  obtained  from  the  prediction  covariance  and  are  related  to  the 

geographic  gate. 


Figure  4  shows  true  target  position  at  sampling  times,  target  measurements,  and 
tracker  corrected  estimates  (zk\k).     In  this  figure,  the  error  ellipses  represent  the  95% 
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confidence  region  around  the  corrected  estimate  and  are  obtained  from  P*]*  .  The  longer 
Lines  of  Bearing  correspond  to  target  2.  The  tracker  is  performing  reasonably  well  with 
both  targets.  Estimate  error  increases  during  target  1  's  turn  but  quickly  diminishes. 
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Figure  4.  Simulation  1  True  Target  Position,  Target  Measurements,  and  Corrected 

Estimates  at  Sample  Times. 


Figure  5  contains  the  plot  of  Figure  4  for  target  1  only.  Likewise,  Figure  6  is  the 
target  2  portion  of  figure  4. 
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Target  1  Motion,  Measurements,  and  Estimates  at  Sample  Times 
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5.  Simulation  1  True  Target  Position,  Target  Measurements,  and  Corrected 
Estimates  at  Sample  Times  for  Target  1. 
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Target  2  Motion,  Measurements,  and  Estimates  at  Sample  Times 
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Figure  6.  Simulation  1  True  Target  Position,  Target  Measurements,  and  Corrected 
Estimates  at  Sample  Times  for  Target  2. 
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Table  1  is  a  compilation  of  statistics  from  the  simulation.    The  following  is  an 
explanation  of  the  fields: 

index:  measurement  number 

report  type:  "1"  indicates  a  cartesian  measurement  sensor  report  (Lat/Long).  "0" 
indicates  a  line  of  bearing  measurement  sensor  report. 

No  of  Clutter  Msmts:  The  number  of  false  measurements  generated  at  that 
measurement  time. 

No  of  Valid  Tracks:  The  number  of  tracks  being  processed  by  the  tracker.  It 

includes  the  original  two  targets  plus  new  tracks  that  were  generated  by  the 
track  splitting  algorithm,  less  tracks  that  were  deleted  for  exceeding  score  or 
maximum  time  without  update  thresholds. 
The  following  statistics  are  compiled  for  both  targets: 

No  w/  Valid  Attributes:  The  number  of  measurements  at  the  sampling  time  that 
have  emitter  frequency  and  PRF  parameters  close  enough  to  the  target's 
expected  values  to  pass  the  feature  rejection  loop. 

No  of  Valid  Measurements:  The  number  of  measurements  at  the  sampling  time  to 
pass  the  feature  rejection  loop  and  fall  within  the  range  gate.  If  this  number 
is  "1",  the  track  is  updated  with  this  measurement.  If  this  number  is  greater 
than  one,  the  measurement  closest  to  the  prediction  will  update  the  track, 
and  remaining  measurements  will  be  used  to  generate  new  tracks  (track 
split).  If  this  number  is  "0",  the  track  will  not  be  updated  at  the  sample  time 
and  output  plots  will  contain  the  predicted,  vice  corrected,  estimate. 
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Msmt  error:  Geometric  distance  in  60  nm  (one  degree)  units  between  the  target 

measurement  and  the  true  target  position. 
Estimate  error:  Geometric  distance  in  60  nm  (one  degree)  units  between  the 

corrected  estimate  and  the  true  target  position. 
Estimate  swap:  The  track  number  of  the  candidate  estimate  with  the  lowest 
likelihood  score.  If  this  number  is  the  same  as  the  target  number,  the 
estimate  for  the  target  is  the  best  available  and  no  swapping  will  occur.  If, 
however,  this  track  number  is  different  from  the  target  number,  the  track 
numbers  will  be  swapped  so  that  the  estimate  for  target  one  has  the  lowest 
likelihood  score. 
No  of  Trk  Candidates:  The  number  of  tracks  that  are  candidates  for  this  target. 
These  include  the  original  track  plus  all  tracks  that  were  generated  by  valid 
measurements  (those  passing  the  feature  rejection  loop  and  falling  within  the 
gate),  less  those  tracks  that  were  deleted  for  exceeding  score  or  maximum 
time  without  update  thresholds. 
Examination  of  the  highlighted  rows  in  Table  1  reveals  the  effectiveness  of  the 
feature  rejection  loop  when  target  and  clutter  frequency  and  PRF  means  are  well  separated 
relative  to  the  sensor's  standard  deviations.   Despite  a  large  number  of  measurements  not 
originating  from  the  target,  few  have  frequency  and  PRF  measurements  close  enough  to  the 
target  expected  values  to  be  considered  valid.  In  the  case  of  the  simulation  results  presented 
here,  typically  only  one  measurement  successfully  passes  this  loop  per  target,  presumably 
originating  from  the  target  being  processed. 
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index 

report  type 

no  of  clutter  msmts 

no  of  valid  tracks 

target  1 
no  of  valid  attributes 
no  of  valid  msmts 
msmt  error 
estimate  error 
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no  of  trk  candidates 
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msmt  error 
estimate  error 
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0 
2 
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1 
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5 
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1 

82 
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7 

1 

42 

2 
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1111111 
0.1065  0.0537  0.1439  0.1056  NaN  0.0908  0.0176 
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1111111 
1111111 

1  1  1  11  1  1 
1111111 

0.0283  0.131 1  0.0162  0.0448  NaN      0.0324  0.0693 
0.0529  0.0977  0.0399  0.0385  0.1040  0.0339  0.0441 

2  2  2  2  2  2  2 
1111111 


8 

1 

48 

2 

1 
1 
0.0490 
0.0329 
1 
1 

1 

0.0874 
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2 
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report  type 

no  of  clutter  msmts 

no  of  valid  tracks 

target  1 
no  of  valid  attributes 
no  of  valid  msmts 
msmt  error 
estimate  error 
estimate  swap 
no  of  trk  candidates 

target  2 
no  of  valid  attributes 
no  of  valid  msmts 
msmt  error 
estimate  error 
estimate  swap 
no  of  trk  candidates 
Table  1 


9 

1 

39 
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0 

39 
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1 
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1 

40 
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2  2  2  2  2  2 
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.  Statistics  for  Simulation  1:  Feature  Rejection  Test 
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1 
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1 
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1 

1 
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0.0937 
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0.0995 

1 

1 

1 

1 

1 

1 

1 

1 
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0.1292 
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0.1378 

2 

2 

1 

1 
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C.         SIMULATION  2:  GATING  TEST 

Another  case  to  consider  is  when  the  two  targets  have  similar  frequency  and  PRF 
means.  This  will  demonstrate  the  effectiveness  of  the  geographic  gating  in  selecting  the 
correct  candidate  measurement.  This  simulation  is  conducted  with  the  same  parameters  as 
the  previous  simulation,  with  the  exception  of  the  target  feature  means: 

f^anl   1.406  GHz 

prfmean-    4.4  ppS 

for  both  targets.  The  false  measurement  feature  means  remain  unchanged,  ensuring  that  the 
majority  of  these  measurements  will  be  rejected. 

Figure  7  shows  true  target  position  at  sampling  times,  target  measurements,  and 
tracker  corrected  estimates  (z^).  As  with  figure  4,  the  error  ellipses  represent  the  95% 
confidence  region  around  the  corrected  estimate  and  are  obtained  from  P^  •  Plots  of  true 
target  track  and  clutter  measurements  have  been  omitted  since  they  are  similar  to  Figures  2 
and  3.  As  expected,  the  tracker  is  performing  reasonably  well  with  both  targets  and 
producing  results  that  are  similar  to  the  case  when  features  are  well  separated. 

Analysis  of  the  simulation  results  contained  in  Table  2  shows  that,  generally,  the 
gating  is  effectively  selecting  the  correct  measurement.  At  each  sample  time,  there  are  two 
measurements  that  have  passed  through  the  feature  rejection  loop,  presumably 
corresponding  to  targets  1  and  2.  Gating  usually  selects  the  correct  measurement  for  use  in 
updating  the  estimate,  preventing  a  track  split.  Occasionally  however,  more  than  one 
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Figure  7.  Simulation  2  True  Target  Position,  Target  Measurements,  and  Corrected 

Estimates  at  Sample  Times. 


measurement  is  in  the  geographic  gate  and  the  track  is  split.  These  instances  are  limited  to 
the  Line  of  Bearing  measurements  (indices  5  and  10),  and  should  be  expected  for  the 
geometry  of  the  simulation.  Referring  to  Figure  3,  the  line  of  bearing  measurements  to  both 
targets  fall  within  the  ellipses  for  both  predicted  estimates,  indicating  a  track  split  situation. 
The  track  is  split  by  updating  the  prediction  with  both  line  of  bearing  measurements,  one  for 
each  track.  The  generated  tracks  have  state  estimates  that  are  very  close  to  the  updated 
track,  making  them  nearly  indistinguishable  in  figure  7. 
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index 

1 

2 

3 

4 

5 

6 

7 

8 

report  type 

1 

1 

1 

1 

0 

1 

1 

1 

no  of  clutter  msmts 

0 

0 

0 

0 

0 

74 

41 

41 

no  of  valid  tracks 

2 

2 

2 

2 

4 

4 

4 

4 

target  1 

no  of  valid  attributes 

2 

2 

2 

2 

2 

2 

2 

2 

no  of  valid  msmts 

1 

1 

1 

1 

2 

1 

1 

1 

msmt  error 

0.1584  0.0286  0.0336  0.1654  NaN 

( 

D.  1298  0.0353 

0.1480 

estimate  error 

0.1 174  0.1008  0.0581  0.1 168  0.2270  0.1458  0.0280 

0.1222 

estimate  swap 

1 

1 

1 

1 

1 

1 

1 

1 

no  of  trk  candidates 

1 

1 

1 

1 

2 

2 

2 

2 

target  2 

no  of  valid  attributes 

2 

2 

2 

2 

2 

2 

2 

2 

no  of  valid  msmts 

1 

1 

1 

1 

2 

1 

1 

1 

msmt  error 

0.0638  0.1803  0.0583  0.1018  NaN 

( 

D.0745  0.0354 

0.0529 

estimate  error 

0.0341  0.1815  0.0528  0.0908  0.1547  0.0751  0.0189 

0.0485 

estimate  swap 

2 

2 

2 

2 

2 

2 

2 

2 

no  of  trk  candidates 

1 

1 

1 

1 

2 

2 

2 

2 

index 

9 

10 

11 

12 

13 

14 

15 

16 

report  type 

1 

0 

1 

1 

1 

1 

1 

1 

no  of  clutter  msmts 

42 

42 

74 

115 

45 

43 

40 

47 

no  of  valid  tracks 

4 

8 

8 

8 

8 

8 

8 

8 

target  1 

no  of  valid  attributes 

2 

2 

2 

2 

2 

2 

2 

2 

no  of  valid  msmts 

1 

2 

0 

1 

1 

1 

1 

1 

msmt  error 

0.0936  NaN 

0.1133  0.1279  0.0443  0.0564  0.1033 

0.0692 

estimate  error 

0.1283  0.3668  0.4907  0.1643  0.0224  0.0620  0.0620 

0.0634 

estimate  swap 

1 

1 

1 

1 

1 

1 

5 

1 

no  of  trk  candidates 

2 

4 

4 

4 

4 

4 

4 

4 

target  2 

no  of  valid  attributes 

2 

2 

2 

2 

2 

2 

2 

2 

no  of  valid  msmts 

1 

2 

1 

1 

1 

1 

1 

1 

msmt  error 

0.1269  NaN 

0.1548  0.1568  0.0887  0.0255  0.0426 

0.2025 

estimate  error 

0.0793  0.1603  0.1 153  0.1263  0.041 1  0.0387  0.0301 

0.1650 

estimate  swap 

2 

2 

2 

2 

2 

2 

2 

2 

no  of  trk  candidates 

2 

4 

4 

4 

4 

4 

4 

4 

Table  2.  Statistics  for  Simulation  2:  Gating  Test 
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D.        SIMULATION  3:  SCORING  TEST 

The  tracker's  ability  to  select  the  best  track  candidate  is  examined  in  a  simulation 
where  the  clutter's  feature  means  are  relatively  close  (within  3  standard  deviations)  of  a 
target's  means: 

Mean  Clutter  Features: 

f clutter-     1.4  GHZ 
prf clutter'-    3.8  ppS 

Mean  Target  1  Features: 
Uanl   1.403  GHz 

Ptf mean'.    4.1  ppS 

The  tracker  will  have  several  false  measurements  successfully  pass  the  feature  rejection 
loop  for  target  1.  Of  these  measurements,  it  is  expected  that  some  will  fall  within  the 
geographic  gate  causing  the  track  to  be  split.  With  successive  measurements,  track  scoring 
is  used  to  select  the  estimate  that  most  closely  conforms  to  the  motion  model  and  eliminate 
unlikely  tracks.  The  number  of  measurements  for  this  simulation  is  increased  from  16  to  20 
to  better  show  the  track  management  capabilities.  Mean  features  for  target  2  are  as 
specified  in  Simulation  1  and  are  expected  to  produce  similar  results. 

Figure  8  shows  true  target  position  at  sampling  times,  target  measurements,  and 
tracker  corrected  estimates  (z^-).  All  tracks  are  shown  here,  including  those  generated  as  a 
result  of  splitting  the  original  tracks.    However,  error  ellipses  are  plotted  only  around  the 
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estimates  most  likely  to  be  targets  1  and  2  (i.e.,  those  two  candidate  tracks  with  the  lowest 
track  scores).  It  is  apparent  that  several  new  tracks  have  been  generated  for  target  1,  and 
that  the  track  with  the  lowest  score  conforms  to  the  true  target  track  well.  Figure  9  shows 
true  target  position  at  sampling  times,  target  measurements,  and  tracker  corrected  estimates 
(zic\k)  for  the  most  likely  track  estimate  only,  illustrating  this  further. 
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Figure  9.  Simulation  3  True  Target  Position,  Target  Measurements,  and  Corrected 
Estimates  at  Sample  Times.  Best  Estimate  Only  For  Target  1. 


An  understanding  of  the  track  management  capability  for  this  tracker  can  be  obtained  by 
examining  Table  3.  Multiple  measurements  possess  valid  features  and  pass  through  the 
feature  rejection  loop  at  many  measurement  times.  Also,  track  splitting  occurs  when  more 
than  one  measurement  falls  within  the  geographic  gate:  samples  9  and  12.  These  new 
tracks  add  to  the  list  of  track  candidates  for  target  1.  The  number  of  track  candidates 
increases  from  1  to  3  at  sample  9  (the  first  split),  then  peaks  at  7  candidates  at  sample  13.  It 
recedes  to  2  candidates  by  the  end  of  the  simulation  as  a  result  of  track  dropping.  By 
comparing  the  track  score  for  all  track  candidates,  the  tracker  identifies  estimates  that  are 


39 


more  likely  to  be  target  1  than  the  currently  held  estimate  on  one  occasion.  At  index  11, 
track  3  is  selected  and  becomes  the  new  track  1.  The  track  shown  in  Figure  7  is  the 
sequence  of  estimates  with  the  lowest  track  score  at  each  sample  time. 
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2 

1  1 
0             0 

2  2 

1 
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2 

0 
0 
2 

1 
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1 

1 
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1              0 

1              0 
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1 
1 

7 
1 

3 
1 

msmt  error 
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0.1495    0.1311 

0.0244  NaN 

0.0640 

0.0719 

estimate  error 
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1 

1 
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1              1 
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1 
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2 

1 
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2 

1 

2 

1 

2 

1 

2 

1 

Table  3. 

Statistics  for  Simulation  3:  Scoring  Test 
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index 
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9 

10 

11 

12 

13 

14 

report  type 

1 

1 

0 

1 
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1 

no  of  clutter  msmts 

45 

43 

43 
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45 

no  of  valid  tracks 

2 

4 

4 

5 

7 

8 

8 

target  1 

no  of  valid  attributes 

3 

5 

5 

7 

5 

5 

7 

no  of  valid  msmts 

1 

3 

1 

0 

2 

1 

1 

msmt  error 

0.0304 

0.0973  NaN 

0.1141 

0.0611 

0.1346 

0.0100 

estimate  error 

0.0269 
1 

0.3519 

1 

0.6189 
1 

0.1039 

0.0890 
1 

0.1273 

1 
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estimate  swap 

3 

4 
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3 
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7 
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1 

1 
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1 

1 

1 
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estimate  error 
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0.1233 
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2 

2 

2 
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2 

2 
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1 

1 
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46 

46 
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47 
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46 
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8 

7 
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3 
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8 

2 

8 

2 

2 
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1 

1 

1 

1 

1 

0 

msmt  error 

0.1348 
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0.0535 
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0.0730 

0.1711 

estimate  error 
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estimate  swap 

1 

1 

1 

1 

1 

1 

no  of  trk  candidates 

7 

6 

3 

2 

2 

2 

target  2 

no  of  valid  attributes 

1 

1 

1 

1 

1 

1 

no  of  valid  msmts 

1 

1 

1 

1 

1 

1 

msmt  error 

0.0684 

0.0740 

0.0996 

0.0962 

0.1436 

0.1052 

estimate  error 

0.0729 

0.0490 

0.0667 

0.1037 

0.0758 

0.0755 

estimate  swap 

2 

2 

2 

2 

2 

2 

no  of  trk  candidates 

1 

1 

1 

1 

1 

1 

Table  3  (Continued).  Statistics  for  Simulation  3:  Scoring  Test 
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Repeated  simulations  with  these  initial  parameters  produce  results  that  vary  in  detail 
but  are  generally  the  same  in  overall  performance.  The  tracker  consistently  obtains  a 
sequence  of  estimates  that  conforms  with  the  true  track  reasonably  well.  Differences 
between  simulations  are  found  in  the  number  and  occurrence  of  track  splitting.  Also,  some 
simulations  have  more  estimate  swaps  than  presented  here,  while  others  have  less.  This  is 
expected  due  to  the  randomness  of  the  clutter  and  target  measurements.  The  track  scoring 
and  comparisons  performed  by  the  tracker  have  demonstrated  considerable  strength  in 
separating  sequences  that  arise  from  uniformly  distributed  clutter  measurements  from  those 
that  arise  from  Gaussian  distributed  target  measurements. 

E.         SIMULATION  4:  COMBINED  TEST 

Finally,  the  case  is  considered  where  false  target  features  are  broadly  distributed 
over  a  range  compared  to  the  target  features.  This  is  the  most  realistic  of  the  scenarios 
presented,  since  it  models  target  tracking  in  an  environment  where  many  other  contacts  are 
operating  with  various  emitters.  A  combination  of  the  tracker's  capabilities  will  be  required 
to  maintain  a  reasonable  track:  feature  rejection,  gating,  and  track  scoring.  Furthermore,  it 
will  be  assumed  that  the  tracker  does  not  know  the  feature  mean  for  the  false  measurements 
and  the  alternate  feature  rejection  technique  discussed  in  Section  B  of  Chapter  HI  will  be 
used.  For  the  purpose  of  generating  false  target  feature  measurements,  a  normal  distribution 
is  used  with  the  following  statistics: 
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Mean  Clutter  Features: 

f clutter-    1.4GHZ 
prf clutter'-    3.8  pps 

Clutter  Standard  Deviation: 

Gfreq.     1.000  GHz 
(Jprf.    3.0  ppS 


Target  features  are  normally  distributed  as  specified  in  Simulation  1 : 
Target  1: 

fmean-.    1.406  GHz 
prfmean-    4.4  ppS 

Target  2: 

fmeanl    1.394  GHz 
prfmean'-    3.8  ppS 

Target  Standard  Deviation: 

(Jfreq-    0.001  GHZ 
(Jprf.    0.1  ppS 
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Figure  10  shows  true  target  position  at  sampling  times,  target  measurements,  and 
tracker  corrected  estimates  (z^).  As  before,  all  tracks  are  shown  here  with  error  ellipses 
plotted  only  around  the  estimates  most  likely  to  be  targets  1  and  2.  The  tracker  is 
performing  reasonably  well  with  both  targets.  Some  new  tracks  have  been  generated 
around  target  1,  but  the  most  likely  state  estimate  conforms  well  with  true  target  position. 


Target  Motion,  Measurements,  and  Estimates  at  Sample  Times 
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Figure  10.  Simulation  4  True  Target  Position,  Target  Measurements,  and  Corrected 

Estimates  at  Sample  Times. 
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Table  4  shows  the  output  statistics  for  this  simulation.  It  is  noted  that  a  few  false 
measurements  pass  through  the  feature  rejection  loop  frequently,  but  track  splitting  occurs 
rarely.  There  is  only  one  track  split  apparent  with  target  1,  and  that  occurs  at  sample  7. 
Target  2  does  not  exhibit  any  track  splitting.  The  number  of  valid  tracks  increases  from  2  to 
3  at  sample  7,  then  from  3  to  4  at  sample  11.  Since  there  is  only  one  valid  measurement  in 
the  gate  for  both  targets  at  this  time,  it  suggests  that  the  track  generated  at  sample  7  was 
split.  The  number  of  valid  tracks  recedes  to  3  at  index  15,  where  it  remains  until  the  end  of 
the  simulation.  It  is  also  noted  that  none  of  the  generated  tracks  ever  achieve  a  lower  score 
than  the  original  two  tracks  and  thus  no  estimate  swapping  occurs. 


index  12            3            4            5  6            7 

report  type  11110  11 

no  of  clutter  msmts  0            0            0            0            0  87          44 

no  of  valid  tracks  2            2            2            2            2  2            3 
target  1 


no  of  valid  attributes 
no  of  valid  msmts 

1 
1 

1 
1 

1 
1 

1 
1 

1 
1 

2 

1 

2 
2 

msmt  error  0.0355   0.1014   0.1158   0.0805  NaN         0.0263   0.0098 

estimate  error  0.0323   0.0898   0.1222   0.0978   0.1060   0.0288   0.0019 


estimate  swap 

1 

1 

1 

1 

1 

1 

1 

no  of  trk  candidates 

1 

1 

1 

1 

1 

1 

2 

target  2 

no  of  valid  attributes  111113           2 

no  of  valid  msmts  111            1111 

msmt  error  0.0421    0.0833   0.0262   0.0819  NaN         0.0804   0.0829 

estimate  error  0.0374   0.0695   0.0372   0.0663   0.1019    0.0723   0.0725 

estimate  swap  2            2            2            2            2            2            2 

no  of  trk  candidates  1111111 
Table  4.  Statistics  for  Simulation  4:  Combined  Test 
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index 

report  type 
no  of  clutter  msmts 
no  of  valid  tracks 
target  1 


8 

1 

42 

3 


9 

1 

42 

3 


10 
0 

42 
3 


11 
1 

81 
4 


12 
1 

48 
4 


13 
1 

40 
4 


14 

1 

43 

4 


no  of  valid  attributes 
no  of  valid  msmts 

2 

1 

2             2 
1             1 

1 
1 

2 

1 

2 
1 

1 

msmt  error 
estimate  error 
estimate  swap 

0.0929 

0.0636 

1 

0.0953  NaN 
0.0945    0.2575 
1              1 

0.0496 

0.0888 

1 

0.0846 

0.0564 

1 

0.0490 

0.0596 

1 

0.1498 
0.0887 

1 

no  of  trk  candidates 

2 

2             2 

3 

3 

3 

3 

target  2 

no  of  valid  attributes 
no  of  valid  msmts 

1 

1 

1             1 
1             1 

1 
1 

1 

1 

1 
1 

] 

msmt  error 
estimate  error 
estimate  swap 

0.1319 

0.0753 

2 

0.0371  NaN 
0.0623   0.0786 
2             2 

0.0363 

0.0276 

2 

0.1444 

0.1021 

2 

0.1112 

0.1122 

2 

0.0291 

0.0518 

2 

no  of  trk  candidates 



1 

1             1 

1 

1 

1 

1 

index 

report  type 

no  of  clutter  msmts 

no  of  valid  tracks 

target  1 

/alid  < 

ralid  i 

error 

Lte  er 

ite  sv 

rkca 

targ< 

ralid  < 

ralid  I 
msmt  error 
estimate  error 
estimate  swap 
no  of  trk  candidates 


15 
1 

46 
3 


16 
1 

44 
3 


17 
1 

42 
3 


18 
1 

50 
3 


19 
1 

45 
3 


20 
1 

45 
3 


no  of  valid  attributes 
no  of  valid  msmts 

1 
1 

1 
1 

1 

1 

2 
1 

2 
1 

1 
1 

msmt  error 
estimate  error 
estimate  swap 

0.0775 

0.0974 

1 

0.0756 
0.0391 

1 

0.0542 
0.0600 

1 

0.0384 
0.0176 

1 

0.1297 

0.0966 

1 

0.1292 

0.1287 

1 

no  of  trk  candidates 

2 

2 

2 

2 

2 

2 

target  2 

no  of  valid  attributes 
no  of  valid  msmts 

1 
1 

1 
1 

1 
1 

2 
1 

2 

1 

1 
1 

0.0436 

0.0329 

2 

1 


0.0706 

0.0600 

2 

1 


0.1265 

0.0749 

2 

1 


0.0961 

0.0536 

2 

1 


0.1294 

0.0859 

2 

1 


0.1203 

0.1296 

2 

1 


Table  4  (Continued).  Statistics  for  Simulation  4:  Combined  Test 
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Just  as  with  the  previous  simulations,  repeated  runs  with  these  initial  parameters 
produce  results  that  vary  in  detail  but  are  consistent  in  overall  performance.  Several  runs 
were  observed  with  no  occurrence  of  track  splitting.  On  the  other  hand,  many  simulations 
produced  satisfactory  results,  but  with  track  splitting  and  estimate  swapping.  The  tracker's 
proven  capabilities  in  splitting,  scoring,  and  swapping  result  in  target  estimates  that  conform 
reasonably  well  with  truth. 


F.         SIMULATIONS  WITH  EXTERNAL  MEASUREMENT  DATA 

Much  of  the  code  contained  in  the  Appendix  is  devoted  to  generating  data  that  the 
tracker  would  receive  from  the  sensor  report.  Further  testing  of  the  algorithm  could  be 
accomplished  by  stripping  off  this  code  and  reading  in  collected  sensor  reports  in  the  form 
of  an  external  data  file.  With  minor  modifications,  the  tracker  could  be  tested  with  position 
reports  on  two  real  tracks  and  several  false  measurements.  Sensor  reports  must  contain  the 
following  information: 


*e 


Time  of  Measurement 

Measurement  Type:  a  flag  for  indicating  Cartesian  (Latitude  and  Longitude)  or 

Line  of  Bearing  processing. 
Position  Measurement 
Error  Ellipse:  The  sensor  covariance  can  be  derived  from  this  as  discussed  in 

Chapter  JJ. 
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Emitter  Frequency 

Frequency  Measurement  Covariance 

Emitter  PRF 

Emitter  PRF  Measurement  Covariance 

Because  the  computer  simulation  in  the  Appendix  is  designed  to  evaluate  the  track 
maintenance  capability  of  a  track  splitting  approach,  care  must  be  taken  in  initializing  the 
tracker.  In  the  Appendix,  the  state  estimates  for  Targets  1  and  2  are  initialized  by  passing 
two  noisy  measurements  of  target  position  to  the  supporting  m-file  kal2init.  The  initial 
state  estimate  and  covariance  are  then  obtained  as  described  in  Chapter  IQ.  A  similar 
approach  could  be  taken  when  using  external  data  by  reading  in  two  measurements  for  each 
target  and  passing  them  to  kal2init.  Another  approach  would  be  to  simply  provide  the 
initial  estimate  as  a  priori  knowledge  and  proceed  directly  to  the  main  loop  in  the  code. 
Similarly,  emitter  frequency  and  PRF  means  are  provided  as  a  priori  knowledge  in  the 
Appendix  code.  When  testing  with  external  data,  these  a  priori  estimates  must  be 
consistent  with  the  data  for  satisfactory  performance. 
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V.         CONCLUSION 

A.         SUMMARY 

This  research  has  developed  a  simulation  to  explore  the  use  of  a  track  splitting 
tracker  in  the  maritime  surveillance  problem.  The  tracker  uses  Kalman  Filtering  to  estimate 
target  state  from  received  measurements,  which  may  be  latitude  and  longitude,  or  line  of 
bearing  reports.  Prior  to  filter  update,  measurements  are  filtered  on  the  basis  of  features  and 
geographic  gates  to  eliminate  unlikely  measurement  candidates.  When  more  than  one 
measurement  falls  within  a  track's  gate,  the  track  is  split  to  allow  subsequent  measurements 
to  help  determine  validity.  Track  scoring  with  a  likelihood  function  serves  as  the  basis  for 
determining  which  track  estimate  originated  from  the  modeled  target's  sequence  of 
measurements.  To  prevent  saturation  of  the  tracker  with  false  tracks,  a  track  dropping 
routine  eliminates  the  least  likely  estimates  by  considering  track  scores  and  time  since  last 
update. 

The  effectiveness  of  measurement  rejection  on  the  basis  of  features  was 
demonstrated  in  simulation  by  assigning  feature  means  to  the  false  measurements  that  were 
well  separated  from  the  true  target  means.  Examination  of  the  simulation  results  showed 
that  few  or  no  clutter  measurements  passed  through  this  rejection  loop,  while  the 
measurements  originating  from  the  target  generally  did.  The  Kalman  Filter  track  estimates 
conformed  well  with  the  true  target  tracks.  To  investigate  the  effectiveness  of  the 
geographic  gating,  both  targets  were  assigned  the  same  feature  means.    Since  the  targets 
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were  well  separated  geographically  but  had  similar  features,  both  target  measurements 
successfully  passed  the  feature  rejection  test  and  were  then  correctly  paired  with  a  track  on 
the  basis  of  the  geographic  gate.  Finally,  the  ability  of  the  tracker  to  select  the  best  track 
estimate  from  the  set  of  candidate  tracks  was  tested  by  assigning  similar  feature  means  to 
the  false  measurements  and  one  of  the  targets.  During  the  simulation,  the  track  was  split 
several  times  as  a  result  of  false  measurements  passing  both  the  feature  rejection  and 
geographic  gate  tests.  The  tracker  produced  an  estimate  of  target  track  that  conformed  well 
with  truth  by  selecting  the  lowest  scored  estimate  at  each  measurement  time.  These 
simulations  demonstrated  several  desirable  features  for  a  maritime  surveillance  tracker. 

B.         FURTHER  RESEARCH 

With  minor  modifications,  the  code  in  the  Appendix  can  be  tested  with  a  real  data 
set.  Initially,  testing  should  resemble  the  simulations  presented  here  with  two  real  targets 
and  several  false  measurements  at  each  measurement  time.  Target  feature  means  must  be 
known  before  hand,  and  target  tracks  must  be  initialized  before  clutter  measurements  are 
introduced.  Such  a  test  would  confirm  the  viability  of  the  track  splitting  tracker  in  a  simple, 
real-world  application. 

The  algorithm  can  be  significantly  enhanced  by  adding  a  more  sophisticated  track 
initialization  routine.  Presently,  the  simulation  tests  the  ability  of  the  track  splitting 
approach  to  maintain  target  tracks  once  obtained.  Measurements  that  do  not  pass  the 
feature  test  and  fall  within  the  geographic  gate  are  discarded.    An  initialization  routine 
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might  establish  tentative  tracks  with  these  measurements,  which  could  then  become  regular 
tracks  when  scoring  or  number  of  update  thresholds  are  met. 

A  more  accurate  geographic  model  can  be  developed  for  this  tracker.  Here, 
distances  are  converted  to  geographic  latitude  and  longitude  with  a  flat  earth  model.  Since 
sensors  may  operate  over  great  distances  in  maritime  surveillance,  a  spherical  earth  model 
should  be  used. 
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APPENDIX 

This  Appendix  contains  the  track-splitting  tracker  simulation  and  supporting  code. 
It  was  implemented  in  MATLAB®  5.1. 

Main  Simulation  Code:  trackSplit.m 

%trackSplit.m 

%     This  file  simulates  a  multi-target,  multi-sensor,  track  splitting  tracker.  Two 
%dimensional  motion  is  simulated  using  a  constant  velocity  discrete  time  target  motion 
%model  for  two  contacts  (Target  1  and  Target  2).  An  instantaneous  turn  may  be 
%specified  by  selecting  a  time  to  turn  and  entering  that  value  for  "turnTime"  below. 
%Targets  also  have  emitter  frequency  and  prf  attributes. 

%      A  measurement  of  target  position  and  attributes  are  obtained  "about  every  90  min". 
%Both  targets  are  measured  at  the  same  instant.  Measurement  noise  is  modeled  as  a  zero 
%mean  Gaussian  Random  Variable  with  variance  as  set  in  the  SENSOR  part  of  the  code. 
%The  position  measurements  may  be  taken  as  Range  and  Bearing  or  Bearing  only. 
%Measurements  are  taken  as  Range  and  Bearing  except  on  measurement  indices  that 
%have  been  specified  for  Bearing  only  measurements.  Also,  an  IF  statement  is  in  place 
%for  use  in  "skipping"  measurements  (missed  measurements). 
%      Clutter  measurements  are  also  generated  with  target  range  and  bearing 
%measurements.  Clutter  measurements  are  modeled  as  uniformly  distributed  within  a 
%region  around  true  track  position.  These  measurements  are  assumed  to  have  frequency 
%and  prf  attributes  that  are  also  modeled  as  Gaussian  random  variables.  These  statistics 
%are  specified  in  the  SENSOR  part  of  the  code.  The  clutter  region  is  related  to  the 
%eigenvalues  of  the  estimate  covariance  matrix,  P(k|k-1),  which  may  extremely  large  for 
%tracks  in  their  infancy.  For  this  reason,  an  IF  statement  controls  when  clutter 
%measurements  begin. 

%      Clutter  rejection  is  accomplished  by  comparing  attributes  of  all  measurements  with 
%the  track's  expected  attributes.  A  weighted  distance  between  clutter  and  target  attribute 
%means  is  calculated  for  each  measurement  and  the  ratio  is  formed.  Measurements  are 
%rejected  as  clutter  when  the  ratio  exceeds  a  threshold,  "attThresh". 
%   From  the  remaining  measurements,  Geometric  Association  is  accomplished  by  a 
%comparison  is  made  by  comparing  expected  target  position  (z(k|k-l)  and  measurement 
%position.  The  measurement  that  produces  a  track  with  the  lowest  log-likelihood  score 
%is  selected  for  filter  update.  New  Tracks  are  generated  for  all  remaining  measurements 
%in  the  validation  region  (Chi-Square  Distribution  95%  Confidence  Ellipse).  A  kalman 
%filter  performs  a  correction  step  using  the  selected  measurement,  resulting  in  z(k|k).  An 
%error  ellipse  is  obtained  from  the  covariance  of  error  matrix  for  plotting. 
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%     Track  Dropping  logic  is  included  to  control  the  number  of  tracks  in  the  system. 

%After  all  tracks  have  been  processed  at  a  given  measurement  time,  track  scores  and 

%last  update  time  are  compared  to  threshold  values  (maxScore  and  maxTime).  Tracks 

%that  exceed  the  threshold  are  dropped. 

%  Output  plots  are: 

%  True  target  motion 

%  True  target  motion  and  measurements 

%  True  target  motion,  measurements,  and  corrected  estimates  (x(k|k))  with 

%  error  ellipse 

%  True  target  motion,  measurements,  and  corrected  estimates  for  target  1 

%  True  target  motion,  measurements,  and  corrected  estimates  for  target  2 


% 

This  simulation  requires  I 

% 

kal2init.m 

% 

elipa.m 

% 

pole2cart.m 

% 

ss 

% 

get2d.m 

% 

% 

getLob.m 

% 

getClut.m 

% 

getClutFeature.m 

ig  additional  m-files: 


{initializes  the  filter  from  2  measurements} 
{ obtains  error  ellipse  from  P  matrix } 
{converts  range  and  bearing  measurements  to 
cartesian  coords } 

{ generates  cartesian  position  reports 
(measurements)} 

{ generates  LOB  position  reports } 

{ generates  "clutter"  cartesian  position  reports } 

{generates  "clutter"  features  reports} 


clear 

%%%%%%%%%%%%%%%Simulation  and  Sensor  Initialization%%%%%%%% 

%Define  number  of  measurements  to  be  taken: 


nmeas  =  16; 
clutDens  =  5; 


%n  umber  of  measurements 

%defines  clutter  density:  5  tgt/3600  sqr  nm 


%sensor 

slp  =  [-113.00;18.00]; 
sigmalr  =  5/60; 
sigmalb  =  .5*pi/180; 


%sensor  location  in  lon;lat  deg 

%2.5  nm  of  uncertainty  in  range 
%l/2  degrees  uncertainty  in  brg 


slv  =  diag([sigmalrA2;sigmalbA2]);  %covariance  matrix  for  uncorrelated 

%range  and  brg  msmts 
sigmaF  =  le-3;  %  100  Hz  uncertainty  in  f 

sigmaPrf  =  0. 1 ;  %0. 1  pps  uncertainty  in  prf 


fClut=  1.400; 
prfClut  =  3.8; 


%clutter  feature  statistical  means 
%Ghz  and  pps 
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sigmaFClut  =  1 ;  %clutter  feature  statistical  STD 

sigmaPrfClut  =  3; 

sigC  =  diag([l/sigmaFClut;  1/sigmaPrfClut]);  %parameter  "packaging"  for  use  in 

sigT  =  diag([l/sigmaF;  1/sigmaPrf]);  %clutter  rejection  loop 

attClut  =  [fClut;prfClut]; 

c  =  sqrt(6);  %error  ellipsoids  for  95%  confidence 

gamma  =  6;  %parameter  for  95%  threshold  of  chi  squared  dist 

attThresh  =  10;  %threshold  for  attribute  discrimination 

maxScore  =  48;  %Track  dropping  parameters 
maxTime  =  270; 

%%%%%%%%%%%%%%%%TargetInitialization%%%%%%%%%%%%%%% 
%Initialize  the  Target  1  state  vector  of  positions  and  velocities 

hdg  =  145*pi/180;  %Target  heading 

speed  =  13/3600;  %Target  speed,  kts 

xi  =  [-1 10.00;  %Initial  x  posit  (deg  Ion) 

speed*  sin(hdg);  %Initial  x  speed  (y/s) 

22.00;  %Initial  y  posit  (deg  lat) 

speed*cos(hdg)];  %Initial  y  speed  (y/s) 

turnTime  =  590;  %turn  590  minutes  into  simulation 

%Attribute  Data  for  Target  1 

fTrue  =  1 .406;  %true  emitter  Freq  (GHz) 

prfTrue  =  4.4;  %true  emitter  prf  (pps) 

attTrack(:,l)  =  [fTrue;prfTrue];  %"packaging" 

%Initialize  the  Target  2  state  vector  of  positions  and  velocities 

hdg2  =  120*pi/180;  %Target  heading 

speed2  =  15/3600;  %Target  speed,  kts 

xi2  =  [-109.00;  %Initial  x  posit  (deg  Ion) 

speed2*sin(hdg2);  %Initial  x  speed  (y/s) 

22.00;  %Initial  y  posit  (deg  lat) 

speed2*cos(hdg2)];  %Initial  y  speed  (y/s) 

turnTime2  =  le6;  %no  turn  for  tgt  2 

%Attribute  Data  for  Target  2 

fTrue2  =  1 .394;  %true  emitter  Freq  (GHz) 

prfTrue2  =  3.2;  %true  emitter  prf  (pps) 

attTrack(:  ,2)  =  [fTrue2;prfTrue2] ;  %"packaging" 


%%%%%%%%%%%%%%Tracker  Parameter  Initialization  %%%%%%%%%%%%% 
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%Defme  the  H  (measurement)  matrix  for  cartesian  measurement 
H=[l  0  0  0; 

0     0  1  0]; 


%Define  Target  1  covariance  of  plant  noise  matrix,  Q 

qsqr  =  le-8;  ^adjustable  to  account  for 

%plant  model  deviations 
del  =  90;  %mean  sample  period 


Q  =  qsqr*  [del  A3/3 

delA2/2 

0          0; 

delA2/2 

del 

0          0; 

0 

0 

delA3/3 

delA2/2; 

0 

0 

delA2/2 

del]; 

%Define  Target  2  covariance  of  plant  noise  matrix,  Q 
qsqr2  =le-8; 
del  =  90; 


Q2  =  qsqr2*[delA3/3    delA2/2 

0          0; 

delA2/2             del 

0          0; 

0                       0 

delA3/3 

delA2/2; 

0                       0 

delA2/2 

del]; 

%Build  the  measurement  time  vector  (msmts  taken  "about"  every  90  min) 

rbMeas  =  90*ones(nmeas-3,l)  +  4*randn(nmeas-3,l);  %range/brg  measmt 

%intervals  in  minutes 

lobl  =  rbMeas(5)-30;  %get  a  line  of  brg  msmt  30  min 

%before  4th  r/b  msmt 

lob2  =  rbMeas(9)-45;  %get  a  line  of  brg  msmt  80  min 

%before  9th   r/b  msmt 

hit  =  [0;  rbMeas([l:4]);  lobl;rbMeas([5:8]);lob2;rbMeas([9:(nmeas-3)])]; 

%"hit"  is  the  vector  of  sample 

%  intervals 
hit  =  round(hit);  %measmt  intervals  in  integer  mins 

%%%%%%%%%%%%%Trackmitialization%%%%%%%%%%%%%%%%%% 

%initialize  filter  1  with  2  measurements 

deltalnit  =  (90  +  4*randn(l));  %sample  interval  for  the 
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%two  msmts  used  in  init 

[xhat(:,l)  phat  xi]  =  kal2init(xi,  H,  Q,  sigmalr,  sigmalb,  deltalnit,  sip); 

phatPack(:,l)  =  reshape(phat,16,l);  %"phatPack"  packages  the 

%covariance  matrix 
%into  a  column  vector 

%initialize  filter  2  with  2  measurements 

[xhat(:,2)  phat2  xi2]  =  kal2init(xi2,  H,  Q2,  sigmalr,  sigmalb,  deltalnit,  sip); 

phatPack(:,2)  =  reshape(phat2,16,l); 


%Initialization  for  tracker 

zout  =  []; 

posout  =  []; 

error  =  []; 

time  =  0; 

zest  =  NaN*ones(4,l); 

x  =  xi; 

ind=l 

xb  =  []; 

yb  =  []; 

xestEll  =  []; 
yestEll  =  []; 
xPredEll  =  []; 
yPredEll  =  []; 
estError=  []; 

y  =  []; 

clutOut  =  []; 
yClutOut  =  []; 

lam  =  zeros(2); 

nTracks=2; 
valTrack  =  [12]; 
trkCand  =  [12]; 


and  Output  Matrices  for  track  1 

%Cartesian  msmts 

%true  cartesian  positions 

%cartesian  measurement  error  (from  truth) 

%sample  times 

%position  estimate  (H*xhat(k|k)) 

%state  initialization 

%sample  index 

%x-coord  for  LOB  plotting 

%y-coord  for  LOB  plotting 

%x-coord  for  corrected  estimate 

%(zest)  95%  ellipses 

%y-coord  for  corrected  estimate 

%(zest)  95%  ellipses 

%x-coord  for  predicted  estimate 

%(H*xPred)  95%  ellipses 

%y-coord  for  predicted  estimate 

%(H*xPred)  95%  ellipses 

%estimate  error  (from  truth) 

%feature  measurement  ([freq;prf]) 

%clutter  cartesian  positions 

%clutter  freature  measurements 

%track  score  and  last  update  time 

%[score;updTime] 

%initial  number  of  tracks 

%initial  vector  of  valid  track  numbers 

%initial  vector  of  track  candidates 


%filter  2  initialization 
zout2  =  []; 
posout2  =  []; 
error2  =  []; 
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time  =  0; 
zest2  =  []; 
x2  =  xi2; 
xe2  =  []; 
ye2  =  []; 
xb2  =  []; 
yb2  =  []; 
xestE112  =  []; 
yestE112  =  []; 
xPredE112  =  []; 
yPredE112  =  []; 
estError2  =  []; 

y2  =  []; 

clutOut2  =  []; 
yClutOut2  =  []; 


%%%%%%%%%%%%%%trackersimulation%%%%%%%%%%%%%%%%%%%% 
%Loop  through  the  simulation,  generating  target  motion  between 
%sample  times  and  taking  measurements  at  hit  times 


for  ii  =  l:length(hit), 
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if  (nTracks  >  50)  break,end; 


%stop  sim  if  capacity  exceeds  50  tracks 


delta  =  hit(ii); 


%time  interval  since  last  measurement: 


%  Define  the  F  Matrix  (Transition  Matrix)  for  discrete  time 
%  target  motion  with  constant  velocity 
F=       [1  delta     0  0; 

0  10  0; 

0  0  1  delta; 

0  0  0  1]; 


Q  =  qsqr*[deltaA3/3 

deltaA2/2 

0           0; 

deltaA2/2 

delta 

0           0; 

0 

0 

deltaA3/3 

deltaA2/2; 

0 

0 

deltaA2/2 

delta]; 

%filter  prediction 
for  trkNo  =  1  mTracks 

phat  =  reshape(phatPack(:, trkNo), 4,4); 

xpred(:,trkNo)=F*xhat(:,trkNo); 
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Ppred=F*phat*F'  +  Q; 
PpredPack(:,trkNo)  =  reshape(Ppred,16,l); 
end 

%obtain  3-sigma  ellipses  around  estimate  for  display  purposes 
if  (ii>5) 

Ppred=reshape(PpredPack(: ,  1  ),4,4); 

K  =  [Ppred(l,l)Ppred(l,3); 

Ppred(3,l)Ppred(3,3)]; 

[xell,  yell,  smaj,  smin,  theta]  =  elipa(K,  c,  xpred(l,l),  xpred(3,l)); 

xPredEll  =  [xPredEll  xell']; 

yPredEll  =  [yPredEll  yell']; 

Ppred=reshape(PpredPack(:  ,2),4,4); 
K  =  [Ppred(l,l)Ppred(l,3); 

Ppred(3,l)Ppred(3,3)]; 
[xell,  yell,  smaj,  smin,  theta]  =  elipa(K,  c,  xpred(l,2),  xpred(3,2)); 
xPredE112  =  [xPredE112  xell']; 
yPredE112  =  [yPredE112  yell']; 
end 


%target  1  motion  step:  x(k+l)  =  F*x(k) 

%target  changes  course  to  180  at  sample  time  before  turnTime 

if  (time(length(time))  >  turnTime) 

x(2)  =  0; 

x(4)  =  -speed; 
end 
x  =  F*x; 

%target  2  motion  step:  x(k+l)  =  F*x(k) 

%target  changes  course  to  180  at  sample  time  before  turnTime 

if  (time(length(time))  >  turnTime2) 

x2(2)  =  0; 

x2(4)  =  -speed2; 
end 
x2  =  F*x2; 


%take  measmts  at  hit  times: 

%take  measurement  from  current  target  state  vector  and  append 

%to  the  true  target  position  output  matrix 
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time  =  [time  time(length(time))+hit(ii)]; 

ztrue  =  H*x; 

posout  =  [posout  ztrue]; 


ztrue2  =  H*x2; 

posout2  =  [posout2      ztrue2]; 
%*****************End  of  i^g^  Motion  /  Filter  Prediction******** 


%obtain  noisy  msmt  from  sensor  1  with  associated  covariance  matrix 

%in  cartesian  coords. 

%The  "if  statement  is  used  to  simulate  randomly  missed  measurements 

if(rand(l)<l), 

sigR  =  sigmalr; 

sigB  =  sigmalb; 

sigma  =  [sigR;  sigB];  %packaging 

%find  difference  in  positions  between  sensor  and  tgt 
diff  =  ztrue-  sip; 
diff2  =  ztrue2  -  sip; 


%target  1  true  bearing  and  range 
r  =  sqrt(diff(l)A2  +  diff(2)A2); 
brg  =  atan2(diff(l),diff(2)); 

%target  2  true  bearing  and  range 
r2  =  sqrt(diff2(l)A2  +  diff2(2)A2); 
brg2  =  atan2(diff2(l),diff2(2)); 

%target  1  attribute  measurements 
fMeas  =  fTrue  +  randn(l)*sigmaF; 
prfMeas  =  prfTrue  +  randn(l)*sigmaPrf; 

y  =  [y  [fMeas ;prfMeas]]; 

%target  2  attribute  measurements 
fMeas2  =  fTrue2  +  randn(l)*sigmaF; 
prfMeas2  =  prfTrue2  +  randn(l)*sigmaPrf; 

y2  =  [y2  [fMeas2;prfMeas2]]; 
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repTyp2d  =  ~((ii=5)|(ii=10));  %specify  msmts  to  be  LOB 

%here,  5th  and  10th  msmt 
%will  be  LOB, 
%all  others  cartesian 

if  repTyp2d,         %bearing  and  range  msmt 

[z(:,l),  R]  =  get2d(r,  brg,  sigma,  sip);  %tgt  1 

%returns  noisy  msmt  and 
%msmt  covariance 

RPack(:,l)  =  reshape(R,4,l); 

ztilde  =  ztrue  -  z(:,l); 

disterror  =  sqrt(ztilde'*ztilde);  %msmt  error 

[z(:,2),  R2]  =  get2d(r2,  brg2,  sigma,  sip);        %tgt  2 

%returns  noisy  msmt  and 
%msmt  covariance 

RPack(:,2)  =  reshape(R2,4,l); 

ztilde2  =  ztrue2  -  z(:,2); 

disterror2  =  sqrt(ztilde2'*ztilde2);%msmt  error 

%add  the  measurement  to  the  measurement  output  matrix 
zout  =  [zout      z(:,l)]; 
zout2  =  [zout2  z(:,2)]; 


%This  block  will  produce  clutter  posit  at  attribute  msmts 
if(ii>5) 

Ppred=reshape(PpredPack(:,  1  ),4,4); 

Sl=H*Ppred*H'  +R; 

[clut,  nClutterl]  =  getClut(Sl,  ztrue,  clutDens,  gamma); 

%obtains  clutter  posit  msmts 
%uniformly  dist  around  tgt  1 

yClut  =  getClutFeature(nClutterl,  fClut,  prfClut,  sigmaFClut, 
sigmaPrfClut);  %obtains  clutter  feature  msmts 

clutOut  =  [clutOut  clut]; 
yClutOut  =  [yClutOut  yClut]; 

Ppred=reshape(PpredPack(:  ,2),4,4); 

S2  =  H*Ppred*H'  +  R2; 

[clut2,  nClutter2]  =  getClut(S2,  ztrue2,  clutDens,  gamma); 
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yClut2  =  getClutFeature(nClutter2,  fClut,  prfClut,  sigmaFClut, 
sigmaPrfClut); 

clutOut2  =  [clutOut2  clut2]; 
yClutOut2  =  [yClutOut2  yClut2]; 

else  %ii<5  ==>  do  not  generate  clutter  msmts 

nClutterl  =  0; 

nClutter2  =  0; 

clut=[]; 

clut2  =  []; 

yClut=[]; 

yClut2  =  []; 
end  %clutter  msmts  for  2d  case 

else  %lob  msmt 

[z(:,l),  xLob,  yLob]  =  getLob(r,  brg,  sigB,  sip); 

%tgt  1  noisy  bearing, 
%and  coordinates  for 
%plotting  LOB 

[z(:,2),  xLob2,  yLob2]  =  getLob(r2,  brg2,  sigB,  sip);  %tgt  2 


%add  the  LOB  measurement  to  the  output  LOB  measurement  matrix 

xb  =  [xb  xLob']; 

yb  =  [ybyLob']; 

xb2  =  [xb2  xLob2'];  yb2  =  [yb2  yLob2']; 

ztilde  =  nan;  %msmt  distance  errors  not  computed 

%  for  LOB 
disterror  =  nan; 
ztilde2  =  nan; 
disterror2  =  nan; 
theta  =  brg; 

end       %measurement  generation 

%add  the  measurement  distance  error  to  measurement  distance  error  matrix 

error  =  [error     disterror]; 

error2  =  [error2  disterror2]; 
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t^******************^,^^  processing  /  Filter  Correction* ************** *** 
%Target  Clutter  Rejection:  Select  correct  measurement  from  vector  of 
%measurements  by  examining  attributes  nClutter  =  nClutterl+nClutter2; 
posVec  =  [z(:,l);  z(:,2);  clut([l:2*nClutterl])';clut2([l:2*nClutter2])']; 

%  vector  of  msmt  positions 

%note  that  the  correct  msmt 

%for  target  1  is  the  1st 

%msmt  in  posVec.  Target 

%2  is  second.  All  others 

%are  clutter 

attVec  =  [fMeas;  prfMeas;  fMeas2;  prfMeas2;  reshape(yClut,2*nClutterl,l); 
reshape(yClut2,2*nClutter2, 1)];  %  vector  of  msmt  attributes 

for  trkNo  =  valTrack  %loops  through  all  tracks 

nVal  =  0;  %initialize  valid  msmt  counter 

Ppred  =  reshape(PpredPack(:,trkNo),4,4); 

%extract  Ppred  for  track  being 
%processed  from  packed  matrix 
forl  =  l:(nClutter+2), 

art  =  [attVec(2*l-l);attVec(2*l)]; 

dT  =  (att  -  attTrackO, trkNo))' *sigT*(att  -  attTrack(:, trkNo)); 
dC  =  (att  -  attClut)'*sigC*(att  -  attClut); 
dTdC(l)  =  dT/dC; 
end  %this  loop  compares  the  distance  of 

%a  msmts  features  to  the  tgt  and 
%clutter  means.  dTdC  is  a  ratio  of 
%these  distances. 

valAtt  =  find(dTdC  <  gamma);  %valAtt  holds  the  index  of  all  msmts 

%that  fall  below  the  threshold. 
%These  msmts  have  attributes  that 
%are  valid  for  this  tgt 
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dTdC=[] ;  %reset  for  next 

%Geometric  Track  Association  and  Filter  Update:  associate  observation  that 

%gives  lowest  score  for  this  track. 

%Track  scoring  determined  by  log-likelihood  function 

%!!!!!!!!!!2  cases!!!!!!!!!!!!!!!!!!! 
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%CASE  1:  2  Dimensional  reports  using  standard  filter 
%CASE  2:  LOB  report  using  EKF 

xhatSplit(:,l)=[lam(l,trkNo)+6;  xpred(:,trkNo);  attTrack(2*trkNo-l); 
attTrack(2*trkNo)]; 

%xhatSplit  is  a  holding  vector  that  will  contain 
%the  score,  state  vector,  and  attributes  that  result 
%from  all  measurement  associations  with  the  track 
%being  processed.  Here,  it  is  initialized  with  the 
%tracks  score,  predicted  state,  and  attributes. 

phatSplitPack(:,l)  =  [lam(l,trkNo)+6;  PpredPack(:,trkNo)]; 

%likewise,  phatSplitPack  holds  the  covariances 
%resulting  from  an  association 

if  repTyp2d  %2d  (cartesian)  report 

R  =  reshape(RPack(:,trkNo),2,2);     %extract  sensor  covariance 
S  =  H*Ppred*H'  +  R;  %S  is  the  innovation  covar 

forl=  l:length(valAtt) 

zObs(:,l)  =  [posVec(2*valAtt(l)-l);posVec(2*valAtt(l))]; 

%obtain  a  valid  (passed  the  attribute  test)  msmt 

ztil  =  zObs(:,l)  -  H*xpred(:,trkNo);      %the  innovation 
lamSplit  =  zuT*inv(S)*ztil;  %norm  of  innovation 

%squared,  a  distance 
%measure 

if  (lamSplit  <  gamma)  %if  the  dist  is  inside  gate 

lam(2,trkNo)  =  time(length(time)); 

%indicate  time  of  update  in 
%score  matrix 
nVal  =  nVal  +  1 ;  %indicate  #  of  valid  associat 

G  =  Ppred*H'*inv(S);  %filter  gain 
xhatSplit(:,nVal)  =  [lam(l,trkNo)+lamSplit; 

(xpred(:,trkNo)  +  G*ztil); 
attVec(2*valAtt(l)-l); 
attVec(2*valAtt(l))]; 
%now,  xhatSplit  holds  new  score  and 
%corrected  est 
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phatSplit  =  (eye(size(Ppred))-G*H)*  Ppred 
*(eye(size(Ppred))-G*H)'+G*R*G' ; 
phatSplitPack(:,nVal)  =  [lam(l,trkNo)+lamSplit; 
reshape(phatSplit,  16,1)]; 
end  %valid  msmt  assoc 
end  %2d  report  association 

else       %LOB  geom  association  and  update 

R  =  sigB;  %sensor  covar 

xl  =xpred(l,trkNo)-slp(l); 
x3  =  xpred(3,trkNo)-slp(2); 
den  =  xlA2  +  x3A2; 

h  =  [-x3/den     0  xl/den  0]; 

%linearized  msmt  matrix 

zhat  =  atan2(xl,x3);  %msmtest 

S  =  h*Ppred*h'  +  R;  %innovation  covar 

forl  =  l:length(valAtt) 

zObs(:,l)  =  [posVec(2*valAtt(l)-l);  posVec(2*valAtt(l))]; 

%obtain  a  valid  observation 

ztil  =  zObs(l,l)  -  zhat;  %innovation 
lamSplit  =  ztil'*inv(S)*ztil;     %norm  of  innovation 

%squared 

if  (lamSplit  <  4)  %if  inside  gate... 

lam(2,trkNo)  =  time(length(time)); 

%indicate  time  of  update 

%in  score  matrix 
nVal  =  nVal  +  1  ;%indicate  #  of  valid  assoc 

G  =  Ppred*h'*inv(S);  %gain 
xhatSplit(:,nVal)  =  [lam(l,trkNo)+lamSplit; 

(xpred(:,trkNo)  +  G*ztil); 
attVec(2*valAtt(l)-l); 
attVec(2*valAtt(l))]; 
%now,  xhatSplit  holds  new  score 
%and  corrected  est 
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phatSplit  =  (eye(size(Ppred))-G*h)  *  Ppred  * 
(eye(size(Ppred))-G*h)'  +G*R*G'; 

phatSplitPack(:,nVal)  =  [lam(l,trkNo)+ 

lamSplit; 

reshape(phatSplit,  16,1)]; 
end 
end 
end       %LOB  geometric  assoc  and  update 

%%% 

if  trkNo=l  %this  IF  block  obtains  the  no 

valAttl  =  length(valAtt);  %of  msmts  that  pass  the 

%feature  clutter  rej- 
nVall  =  nVal;  %ection  and  geographic  gates 


%for  tracks  1  and  2. 


elseif  trkNo==2 


%These  values  will  be  passed 
%to  the  matrix  STAT 
valAtt2  =  length(valAtt);  %for  examination  in  the 

nVal2  =  nVal;  %command  window  after 

%simulation  is  concluded 


end 


%after  simulation,  this  info  can  be  reviewed  in  the  command 
%window.  For  trk  1,  the  index  "1"  should  be  in  the  valAtt  vector.  If 
%not,  attThresh  may  be  too  high.  Likewise,  index  "2"  should  be  included 
%in  valAtt  for  trkNo  2. 

xhatSplit  =  sortrows(xhatSplit');%sorts  validated  estimates  into 
attSplit  =  xhatSplit(:,[6:7])';     %ascending  order  of  likelihood 

%function.  attSplit  holds  the 

%attribute  msmt 

phatSplitPack  =  sortrows(phatSplitPack');  %sorting  by  likelihood 

%here  keeps  covariances 
%indexed  with  estimates 

xhat(:,trkNo)  =  xhatSplit(l,[2:5])'; 
phatPack(:, trkNo)  =  phatSplitPack(l,[2:17])'; 
lam(l,  trkNo)  =  xhatSplit(l,l); 

%lowest  score  is  chosen  for  tgt  track. 

%Estimate,  Covar,  and  score  are 

%updated. 
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for  1  =  2:nVal  %all  remaining  associations 

d  =  size(xhat);  %initialize  new  tracks. 

xhat(:,d(2)+l)=xhatSplit(l,[2:5])';        %...new  trk  state  est 
phatPack(:,d(2)+l)=phatSplitPack(l,[2:17])';  %...new  trk 

%est  covar 
lam(l,d(2)+l)=xhatSplit(l,l);  %...new  trk  score 

lam(2,d(2)+l)=time(length(time));       %...new  trk  init  time 
attTrack(:,d(2)+l)=attSplit(:,l);  %...new  trk  att 

RPack(:,d(2)+l)  =  RPack(:,trkNo);      %...new  trk  sensor 

%covar 
zest(2*d(2)+l,:)  =  nan;  %...new  track  in  estimate 

zest(2*d(2)+2,:)  =  nan;  %output  matrix 

trkCand(d(2)+l)  =  trkCand(trkNo); 
valTrack  =  [valTrack  d(2)+l]; 

end 

xhatSplit  =  [];  %returns  this  to  a  column  vector 

phatSplitPack  =  []; 
end  %trk  processing  loop 

%%%%Track  Management  Algorithm:  2  functions%%% 

%1)  Drop  Tracks  that  exceed  thresholds  (score  and  last  update). 

%2)  Select  candidate  with  lowest  score  to  be  target  estimate 

%  After  all  tracks  have  been  processed  for  a  sample  time... 

%track  dropping... 

delTrk  =  find((lam(l,:)>maxScore)  |  (time(length(time))-lam(2,:))>maxTime); 

%form  a  vector  of  track  numbers  for  track  scores 

%that  exceed  maxScore  or  update  times  that  exceed 

%maxTime 
ifdelTrk 

forl  =  delTrk 

valTrack(find(valTrack==l))  =  []; 

%deletes  track  from  vector  of 
end  %valid  tracks 

end 

xhat(:  ,delTrk)=nan; 

phatPack(:,delTrk)=nan;  %Not-a-number  indicates  "track  deleted" 

lam(:  ,delTrk)=nan; 

%best  track  selection... 

%this  gets  the  track  number  of  the  candidate  w/  lowest  score 
newl  =  floor(find(lam==min(lam(l,find(trkCand==l))))/2)+l; 
new2  =  floor(find(lam=min(lam(l,find(trkCand==2))))/2)+l; 
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%these  lines  swap  the  state  vector  and  co variance  matrices 

%to  put  the  "new"  (lower  score)  vectors  into  trkNo  1  position. 

oldxhat=xhat(:,l); 

oldphat=phatPack( : ,  1 ) ; 

oldlam=lam(:,l); 

xhat(:,l)  =  xhat(:,newl); 
phatPack(:,l)  =  phatPack(:,newl); 
lam(:,l)  =  lam(:,newl); 

xhat(:,newl)  =  oldxhat; 
phatPack(:,newl)  =  oldphat; 
lam(:,newl)  =  oldlam; 

%... likewise  for  trkNo  2 
oldxhat=xhat(:,2); 
oldphat=phatPack(:  ,2); 
oldlam=lam(:,2); 

xhat(:,2)  =  xhat(:,new2); 
phatPack(:,2)  =  phatPack(:,new2); 
lam(:,2)  =  lam(:,new2); 

xhat(:,new2)  =  oldxhat; 
phatPack(:,new2)  =  oldphat; 
lam(:,new2)  =  oldlam; 


else  %missed  msmt 

error  =  [error  nan]; 
end 

d  =  size(xhat); 

nTracks  =  d(2);  %#  columns  in  xhat  is  #  of  tracks  held 

%add  estimate  to  the  matrix  of  estimates 
for  l=l:nTracks 

ze(:,l)  =  H*xhat(:,l); 
end 

zest  =  [zest  reshape(ze,2*l,l)]; 
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%obtain  target  1  error  ellipse  of  estimate 
phat  =  reshape(phatPack(:,l),4,4); 
K  =  [phat(l,l)phat(l,3); 
phat(3,l)phat(3,3)]; 

[xell,  yell,  smaj,  smin,  theta]  =  elipa(K,  c,  xhat(l,l),  xhat(3,l)); 
xestEll  =  [xestEll  xell']; 
yestEll  =  [yestEll  yell']; 

ze=[zest(l,length(zest(l,:)));  zest(2,length(zest(2,:)))]  ; 
estError  =  [estError  sqrt((ztrue-ze)'*(ztrue-ze))]; 

%obtain  target  2  error  ellipse  of  estimate 
phat2  =  reshape(phatPack(:,2),4,4); 
K  =  [phat2(l,l)phat2(l,3); 
phat2(3,l)phat2(3,3)]; 

[xell,  yell,  smaj,  smin,  theta]  =  elipa(K,  c,  xhat(l,2),  xhat(3,2)); 
xestE112  =  [xestE112  xell']; 
yestE112  =  [yestE112  yell']; 

ze=[zest(3,length(zest(3,:)));  zest(4,length(zest(4,:)))] ; 
estError2  =  [estError2  sqrt((ztrue2-ze)'*(ztrue2-ze))]; 

%compile  statistics  for  output 

stat(l,ii)=ii;  %the  matrix  STAT  will  be  displayed  in  the 

stat(2,ii)=repTyp2d;  %command  window  for  post-simulation 

stat(3,ii)=nClutter;  %performance  analysis 

stat(4,ii  )=length(  valTrack) ; 

stat(5,ii)=nan; 

stat(6,ii)=valAttl; 

stat(7,ii)=nVall; 

stat(8,ii)=disterror; 

stat(9,ii)=estError(length(estError)); 

stat(10,ii)=newl; 

stat(  1 1  ,ii)=length(fmd(trkCand(valTrack)==  1 )); 

stat(12,ii)=nan; 

stat(13,ii)=valAtt2; 

stat(14,ii)=nVal2; 

stat(15,ii)=disterror2; 

stat(  1 6,ii)=estError2(length(estError2)); 

stat(17,ii)=new2; 

stat(18,ii)=length(find(trkCand(valTrack)==2)); 


end  %time  step 

%%%%%%%%%%%%%%%OutputResults%%%%%%%%%%%%%%%% 
time(l)  =  []; 
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%Plot  the  results 

nEst  =  nmeas; 

nmeas  =  max(size(zout)); 

%Target  motion  and  true  position  at  measurement  times 

figure(l) 

axis  equal 

plot(posout(l,:),posout(2,:),'-',posout2(l,:),posout2(2,:),'-',slp(l),slp(2),'X');grid; 

title(True  Target  Motion'); 

xlabel('Deg  Lon');ylabel('Deg  Lat') 

%Target  motion  with  noisy  measurements 

figure(2) 

nbrg  =  min(size(xb)); 

nPred  =  min(size(xPredEll)); 

plot(posout(l,:),posout(2,:),'o',zout(l,:),zout(2,:),'*',xb(:,[l:nbrg]),yb(:,[l:nbrg]), 

xb2(:,[l:nbrg]),yb2(:,[l:nbrg]),clutOut(l,:),clutOut(2,:),'*',slp(l),slp(2),'X'); 

hold  on 

plot(posout2(l,:),posout2(2,:),'o',zout2(l,:),zout2(2,:),'*',clutOut2(l,:),clutOut2(2,:),'*', 
slp(l),slp(2),'X,,xPredEll(:,[l:nPred]),yPredEll(:,[l:nPred]), 
xPredE112(:,[l:nPred]),yPredE112(:,[l:nPred]),slp(l),slp(2),'X');grid; 

axis  equal 

title( 'Target  Motion  and  Measurements  at  Sample  Times'); 

xlabel('Deg  Lon');ylabel('Deg  Lat') 

hold  off 

%Target  motion  with  Filter  Estimates 

figure(3) 

plot(posout(l,:),posout(2,:),'o',zout(l,:),zout(2,:),'*',xb(:,[l:nbrg]),yb(:,[l:nbrg]), 

xb2(:,[l:nbrg]),yb2(:,[l:nbrg]),zest(2*[l:nTracks]-l,:),zest(2*[l:nTracks],:),'x', 
xestEll(:,[l:nEst]),yestEll(:,[l:nEst]),slp(l),slp(2),'X');grid; 

hold  on 

plot(posout2(l,:),posout2(2,:),'o',zout2(l,:),zout2(2,:),'*\ 

xestE112(:,[l:nEst]),yestE112(:,[l:nEst]),slp(l),slp(2),'X');grid; 

axis  equal 

title('Target  Motion,  Measurements,  and  Estimates  at  Sample  Times'); 

xlabel('Deg  Lon');ylabel('Deg  Lat') 

hold  off 

%Target  1  motion  with  Filter  Estimates 

figure(4) 

plot(posout(l,:),posout(2,:),'o',zout(l,:),zout(2,:),'*',xb(:,[l:nbrg]),yb(:,[l:nbrgl), 

zest(l,:),zest(2,:),'x',xestEll(:,[l:nEst]),yestEll(:,[l:nEst]),slp(l),slp(2),'X');grid; 
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axis  equal 

title( 'Target  1  Motion,  Measurements,  and  Estimates  at  Sample  Times'); 

xlabel('Deg  Lon');ylabel('Deg  Lat') 

%Target  2  motion  with  Filter  Estimates 
figure(5) 

plot(posout2(l,:),posout2(2,:),'o',zout2(l,:),zout2(2,:),'*,,xb2(:,[l:nbrg]),yb2(:,[l:nbrg]),ze 
st(3,:),zest(4,:),'x',xestE112(:,[l:nEst]),yestE112(:,[l:nEst]),slp(l),slp(2),'X');grid; 
axis  equal 

title( 'Target  2  Motion,  Measurements,  and  Estimates  at  Sample  Times'); 
xlabel('Deg  Lon');ylabel('Deg  Lat') 


stat 


kaOinitm 

function  [xhat,  phat,  xiNew]  =  kal2init(xi,  H,  Q,  sigR,  sigB,  delta,  sensPos) 

%initializes  the  kalman  filter  based  on  2  measurements 

% 

%The  2  "true"  positions  are  xi  and  F*xi. 

%These  are  converted  to  noisy  measurements  using  pol2cart.  The  associated 

%cartesian  covariance  matrices  are  also  obtained. 

%From  these  noisy  measurements  and  covariances,  the  initial  state  and  covariance 
%for  the  filter  is  estimated. 

%ESfPUT  parameters: 


% 

xi: 

true  target  state  at  time  zero 

% 

H: 

Measurement  matrix 

% 

Q: 

Plant  Noise 

% 

sigR: 

Sensor  Range  Inaccuracy  (std  dev) 

% 

sigB: 

Sensor  Bearing  Inaccuracy  (std  dev) 

% 

delta: 

Sample  Interval 

% 

sensPos: 

Sensor  Position 

%OUTPUT  parameters 

%         xhat:  Initial  State  Estimate 

%         phat:  Initial  State  Estimate  Covariance 

%         xiNew:  new  true  target  state  (xi  at  time  delta) 

%Supporting  m-files: 
%         pole2cart.m 
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%  Define  the  F  Matrix  (Transition  Matrix)  for  discrete  time 
%  target  motion  with  constant  velocity 
F  =  [l   delta     0  0; 

0   1  0  0; 

0  0  1  delta; 

0  0  0  1]; 


[zl  Rl]  =  pole2cart(H*xi,  sensPos,  sigR,  sigB); 

xiNew  =  F*xi; 

[z2  R2]  =  pole2cart(H*xiNew,  sensPos,  sigR,  sigB); 


z  =  [z2;  zl]; 

A  =  [l 

1/delta 

0 

0 


0  0  0; 

0  -1/deltaO; 

1  0  0; 
1/delta  0  -1/delta]; 


invf      =  inv(F); 

xhat  =  A*z; 

phat  =  A*[R2  zeros(2,2);zeros(2,2)    Rl+H*invf*Q*invf*HT]*A'; 


pole2cart.m 

function  [z,  R,  disterror]  =  pole2cart(ztrue,  sp,  sigmar,  sigmab); 


%obtains  a  noisy  polar  measurement  to  target,  adds  measurement  noise, 
%and  converts  measurement  and  covariance  matrix  to  cartesian 
%inputs: 

true  posit  of  target  in  cartesian 
posit  of  sensor  in  cartesian 
range  msmt  standard  deviation 

bearing  msmt  standard  deviation 


noisy  measurement  in  cartesian 
covariance  matrix  in  cartesianj 
=  sqrt((ztrue  -  z)*(ztrue  -  z)) 


% 

ztrue  =  [ 

st;yt] 

% 

sp 

;  [xs;  ys] 

% 

sigmar  = 

scalar 

% 

sigmab  = 

=  scalar 

%outputs: 

% 

z  =  [xm; 

ym] 

% 

R  =  2X2 

% 

disterror 
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sigma  =  [sigmar;  sigmab] ; 

%find  difference  in  positions  between  sensor  and  tgt 
diff  =  ztrue  -  sp; 

%true  bearing  and  range 
r  =  sqrt(diff(l)A2  +  diff(2)A2); 
brg  =  atan2(diff(l),diff(2)); 
zplr=[r;brg]; 

%add  noise  to  the  measurement 
zmeas  =  zplr  +  randn(size(sigma)).*sigma; 

%convert  the  noisy  measurement  to  cartesian 
z  =  [zmeas(l)*sin(zmeas(2)); 
zmeas(  1  )*cos(zmeas(2))] ; 
z  =  sp  +  z; 

%convert  covariance  matrix  to  cartesian 
sv  =  diag(sigma.A2); 
M  =  [sin(zmeas(2))   -zmeas(l)*cos(zmeas(2)); 
cos(zmeas(2))       zmeas(l)*sin(zmeas(2))]; 

R  =  M*sv*M'; 


%calculate  distance  error 

ztilde  =  ztrue  -  z; 

disterror  =  sqrt(ztilde'*ztilde); 


elipa.m 

function  [xout,  yout,  smaj,  smin,  theta]  =  elipa(PK,  c,  xt,  yt) 

%calculates  error  ellipsoids  given  error  covariance 
%and  estimate  position 

%INPUTS: 

%         PK  is  covariance  matrix 

%         c  is  conficence  region  (c=sqrt(6)  for  95%) 

%         xt  is  ellipse  center  x  coord 

%         yt  is  ellipse  center  y  coord 
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%OUTPUTS: 

%         xout  is  a  vector  containing  x  coord  points  of  the  ellipse 
%         yout  is  a  vector  containing  y  coord  points  of  the  ellipse 
%         smaj  is  the  length  of  the  semi-major  axis 
%         smin  is  the  length  of  the  semi-minor  axis 
%         theta  is  the  orientation  of  the  semi-major  axis 

%  to  plot  the  ellipse,  use:  plot(xout,yout) 

%adapted  from  Stephen  L.  Spehn's  Errellip.m  (15  Nov  89)  by  Mark  Olson 

%get  eigenvalues  (lam)  and  eigenvectors(V) 
[V,lam]  =  eig(PK); 

sigx  =  sqrt(lam(l,l)); 

sigy  =  sqrt(lam(2,2)); 

%parameterized  ellipse 
t  =  0:2*pi/100:2*pi; 
x  =  sigx*c*cos(t); 
y  =  sigy*c*sin(t); 

%translate  to  eigenvectors  space  and  center  at  tgt  posit 
xout  =  x*V(l,l)  +  y*V(l,2)  +  xt; 
yout  =  x*V(2,l)  +  y*V(2,2)  +  yt; 

%report  semimajor  and  semiminor  axes  lengths  and  orientation 

smaj  =  2*c*max([sigx  sigy]); 

smin  =  2*c*min([sigx  sigy]); 

theta  =  atan2(V(2,l),V(l,l))  +  pi/2*(sigy  >  sigx);        %sigy  >  sigx  =>  y  is  smaj 


get2d.m 

function  [z,  R]  =  get2d(r,  brg,  sigma,  sensPos) 

%this  function  generates  a  noisy  cartesian  measurement  from 
%true  bearing  and  range.  Polar  measurement  noise  is  added  to 
%true  polar  position,  then  rotated  into  a  cartesian  coordinate 
%reference.  The  measurement  covariance  is  also  rotated  into  a 
%cartesian  coordinate  reference. 

%INPUTS: 

%  r:  true  range  to  target 
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% 

brg: 

% 

sigma: 

% 

sensPo 

%OUTPUTS: 

% 

z: 

% 

R: 

true  bearing  to  target 

column  vector  of  sensor  std  dev  [sigB;sigR] 

sensor  position  in  cartesian  coordinates 


cartesian  measurement 

cartesian  measurement  covariance 


%obtains  cartesian  measurement  true  range  and  brg  msmt 

zplr=[r;brg]; 

%add  noise  to  the  measurement 

zmeas  =  zplr  +  randn(size(sigma)).*sigma; 

%convert  the  noisy  measurement  to  cartesian  for  plotting 
z  =  [zmeas(l)*sin(zmeas(2)); 
zmeas(  1  )*cos(zmeas(2))] ; 
z  =  sensPos  +  z; 


%obtain  target  1  covariance  of  error  in  cartesian 

M  =  [sin(zmeas(2))      -zmeas(l)*cos(zmeas(2)); 
cos(zmeas(2))  zmeas(l)*sin(zmeas(2))]; 

R  =  M*diag(sigma.A2)*M'; 


getLob.m 

function  [zbrg  ,xLob,  yLob]  =  getLob(r,  brg,  sigB,  sip); 

%this  function  generates  a  noisy  measurement  bearing  and  output 
%vectors  for  plotting  from  true  target  position. 

%INPUTS 

%         r:  true  target  range 

%         brg:  true  target  bearing 

%         sigB:  sensor  bearing  measurement  std  dev 

%         sip:  sensor  position  in  cartesian  coord 

%OUTPUTS 

%         zbrg:  measured  bearing 
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%         xLob:  Line  of  Bearing  x-coords  for  plotting 

%         yLob:  Line  of  Bearing  y-coords  for  plotting 

%to  plot  the  Line  of  Bearing,  use 
%         plot(xLob,yLob) 

%find  Line  of  Bearing  Msmt  to  target 

%add  noise  to  the  measurement 
zbrg  =  brg  +  randn*sigB; 

%convert  the  noisy  measurement  to  cartesian  for  plotting 
xLob  =  linspace(3,(r+l)*sin(brg)); 
yLob  =  xLob/tan(zbrg); 
xLob  =  slp(l)*ones(size(xLob))  +  xLob; 
yLob  =  slp(2)*ones(size(xLob))  +  yLob; 


getCIutm 


function  [clut,  nClutter]  =  getClut(S,  ztrue,  clutDens,  gamma); 

%this  function  generates  uniformly  distributed  clutter  position  measurement 
%in  a  region  around  the  true  target  position.  (Bar  Shalom  and  Li,  1995) 


Innovation  Covariance  S  =  H*P*H'  +  R 
true  cartesian  position  of  target,  [x-coord;y-coord] 
Clutter  Density,  tgts/area  (here,  tgts/3600  sq  nm) 
chi-squared  dist  threshold  corresponding  to  validation  region 
size,  (here,  gamma=sqrt(6)  for  95%  probability  ellipse 


matrix  of  clutter  positions  clut  =  [x-coord  of  all  clutter  msmts; 

y-coord  of  all  clutter  msmts] 
number  of  clutter  msmts  generated 

%obtain  clutter  posit  measmts 

%generate  target  1  clutter  measurements:  uniformly  distributed  in  a 
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%INPUTS: 

% 

S: 

% 

ztrue: 

% 

clutDens: 

% 

gamma: 

% 

%OUTPUTS: 

% 

clut: 

% 

% 

nClutter: 

%window  around  the  true  target  position: 

%determine  the  area  of  the  validation  region  (ref:  YBS  MMT  sec  3.4.6) 
Vk  =  pi*gamma*sqrt(det(S)); 

nClutter  =  round(10*Vk*clutDens+l);%Number  of  false  msmsts  in  a  square 

%around  truth 

xClutter  =  sqrt(10*Vk)*(rand(l, nClutter)  -  .5*ones(l, nClutter))  +  ztrue(l); 
yClutter  =  sqrt(  10*  Vk)*(rand(l, nClutter)  -  .5*ones(l, nClutter))  +  ztrue(2); 

clut  =  [xClutter;  yClutter]; 


getClutFeaturcm 

function  yClut  =  getClutFeature(nClutter,  fClut,  prfClut,  sigmaFClut,  sigmaPrfClut); 

%this  file  obtains  the  feature  (emitter  freq  and  prf)  measurements  for 
%the  false  targets.  Normal  dist  is  assumed. 


the  number  of  clutter  targets 
mean  emitter  frequency  for  clutter  targets 
mean  emitter  prf  for  clutter  targets 
std  dev  of  emitter  freq  for  clutter  targets 
std  dev  of  emitter  prf  for  clutter  targets 


%INPUTS: 

% 

nClutter: 

% 

fClut: 

% 

prfClut: 

% 

sigmaFClut: 

% 

sigmaPrfClut: 

%OUTPUT: 

% 

yClut: 

% 

% 

matrix  of  feature  measurement  column  vectors 
[emitter  freq  for  all  targets; 
emitter  prf  for  all  targets] 


fMeasClut  =  fClut*ones(l,  nClutter)  +  randn(l, nClutter)*  sigmaFClut; 
prfMeasClut  =  prfClut*ones(l, nClutter)  +  randn(l,nClutter)*sigmaPrfClut; 

yClut  =  [fMeasClut;prfMeasClut]; 
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